Modelado y estudio de un actuador SDBD para el control de flujo. Caracterización de la fuerza electrohidrodinámica generada

Equation Chapter 1 Section 1 Proyecto Fin de Carrera Ingeniería Aeroespacial Modelado y estudio de un actuador SDBD para el control de flujo. Caracterización de la fuerza electrohidrodinámica generada. Autor: Miguel Ángel Domínguez Jiménez Tutor: Miguel Pérez-Saborid Sánchez-Pastor Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 2 Proyecto Fin de Carrera Ingeniería Aeroespacial Modelado y estudio de un actuador SDBD para el control de flujo. Caracterización de la fuerza electrohidrodinámica generada. Autor: Miguel Ángel Domínguez Jiménez Tutor: Miguel Pérez-Saborid Sánchez-Pastor Profesor Titular Dpto. de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 4 Proyecto Fin de Carrera: Modelado y estudio de un actuador SDBD para el control de flujo. Caracterización de la fuerza electrohidrodinámica generada. Autor: Miguel Ángel Domínguez Jiménez 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 6 A mi familia, pareja y amigos. 8 Agradecimientos Con este proyecto pongo fin la etapa de grado, una época inolvidable, que ha supuesto un esfuerzo enorme y muchos sacrificios pero que ha merecido la pena a cada instante. Además de la sólida base de conocimientos que he adquirido, termino siendo una persona mucho más madura, culta y racional de lo que era cuando llegué el primer día a la Escuela. Llegados a este punto y viéndolo todo desde otra perspectiva, no es lo duro y lo exijente del grado lo que me llevo, ni ya me pesan las interminables horas sentado en mi escritorio. Cuando pienso en mi paso por la Escuela recuerdo a cada una de las personas que me han acompañado y los momentos felices que he vivido con ellos. En primer lugar, mis padres. A ellos quiero agradecerles todo lo que han tenido que hacer para que jamás me haya faltado nada y también la educación y la mentalidad luchadora que desde pequeño me han inculcado. Mi padre siempre me cuestionaba, siempre me decía: “Yo creo que no tienes capacidad para eso” o “Que poquito vas a durar en la carrera”. Lo hacía queriendo, porque sabe que confio mucho en mí y en mis capacidades, y yo sin darme cuenta trabajaba mucho más, mi padre “picándome” conseguía ser el motor de mis esfuerzos. Me he esforzado en cada instante de este grado por poder decirle “¡Te lo dije!, Y en cuatro años...”. Sin embargo, ahora sólo quiero decirle “Gracias Papá”. Por otro lado, mi madre, ella siempre ha confiado en mí y me lo ha hecho saber, con ella me reía siempre cuando se iba mi padre y decíamos “Luego se va a tener que callar…”. Ella ha sido un gran apoyo durante estos años, y lo seguirá siendo siempre, así que de nuevo: “Gracias Mamá”. Tampoco me olvido de mi hermana, la muy ingenua se ha metido en la misma carrera, no sabe la que le espera… Sin embargo, me alegro mucho por ella, se que va a crecer muchísimo como persona y que supodrá, como para mí, una etapa inolvidable. Gracias Miriam por siempre estar ahí, incluso cuando esto me ha vuelto insoportable. Y si hablo de soportar, hablo de Mariola, ella lleva aguantando mis cabreos, mis cambios de humor, mis meses encerrado en casa y sobretodo mi falta de tiempo durante todos estos años. Ella, sin tener ninguna necesidad, ha acabado sacrificándose tanto como yo. Por ello, para mí se lo merece todo, porque estoy seguro de que sin ella nada sería lo mismo y que el camino se me habría hecho infinitamente más largo. Por último, y no menos importante, quiero agradecer a mis hermanotes: Andrés, Josema y Pablo que apareciesen en mi vida. Ellos han sido el principal motivo por el cual iba a clase con ganas. Andrés fue la primera persona con la que me cruce en el grado, y sin duda alguien que ha influido muchísimo en mí y en el que confío como en nadie. Pablo es Pablo, es muy suyo y muy especial, y eso lo hace un tío increíble, muy rial. Es una de esas personas que dejan huella y de las que no querría distanciarme. Por último, el que llegó más tarde a la hermandad, Josema. Lo que pasa por su mente asusta, tiene una visión de las cosas completamente distinta a la del resto de personas y nos complementa a los otros tres a la perfección. Los momentos que hemos vivido tanto en Escuela, como en la Feria, en Sevilla, en la playa o simplemente en una videollamada han sido los recuerdos más felices y divertidos que me llevo. Sin lugar a dudas, se han convertido en mis mejores amigos, y gracias a cómo entre todos nos hemos respaldado hemos conseguido llegar a este punto de una manera mucho más fácil y agradable. Agradezco también el trabajo y dedicación del profesorado, especialmente a algunos docentes como Gañan, Sergio Esteban o Carrión, con los que pude establecer una relación muy buena y que me ayudaron bastante. Mención especial a Miguel Pérez-Saborid, con él hemos tenemos una muy buena amistad. Incluso solucionado dudas en el despacho, tanto Pablo como yo nos hemos reido muchísimo con él. Y ya que estamos … Saborid, aunque no lo sepas, Pablo y yo te llamamos “El tito Saborid” porque siempre has estado dispuesto a ayudarnos y a pasar el tiempo que fuese necesario para solucionar nuestros problemas, así que gracias Miguel. Miguel Ángel Domínguez Jiménez 16 ÍNDICE DE ILUSTRACIONES Ilustración 1 – A la izquierda estructura de un gas, a la derecha estructura de un plasma. 17 Ilustración 2 – Descarga de plasma, ejemplo de comportamiento colectivo. 18 Ilustración 3 – Nebulosa de Orion, región de gas ionizado, todos los plasmas emiten luz. 18 Ilustración 4 – Flujo inducido por el actuador de plasma. 20 Ilustración 5 – Relación entre el voltaje y la intensidad en diferentes descargas eléctricas en gases. 21 Ilustración 6 – Descarga de Barrera Dieléctrica Filamentosa [28]. 22 Ilustración 7 – Descarga de Barrera Dieléctrica Difusa. 23 Ilustración 8 – Actuadores DBD superficiales, con descarga única y doble. 23 Ilustración 9 – Material dieléctrico sin y con presencia de un campo eléctrico. 24 Ilustración 10 – Arquitectura básica del actuador. 25 Ilustración 11 – Canal tipo DBD de flujo impulsado por plasma. 26 Ilustración 12 – Dominio y condiciones de contorno para resolver la ecuación (14). 29 Ilustración 13 – Dominio y condiciones de contorno para calcular la densidad de carga. 31 Ilustración 14 – Medida 3D de la emisión de luz durante un ciclo completo de corriente AC [17]. 31 Ilustración 15 – Curvas de nivel de la emisión de luz durante un ciclo completo de corriente AC [17]. 32 Ilustración 16 – Distribución semi-Gaussiana para la densidad de carga. 32 Ilustración 17 – Condiciones de contorno del problema fluido. 33 Ilustración 18 – Condiciones de contorno del problema fluido adimensional. 34 Ilustración 19 – Función de corriente. 34 Ilustración 20 – Condiciones de contorno del problema fluido usando la función de corriente. 35 Ilustración 21 – Regiones consideradas en el problema. 36 Ilustración 22 – Mallado del dominio en x. 36 Ilustración 23 - Aproximación polinómica con nodos equiespaciados y nodos de Chebyshev. 37 Ilustración 24 – Mallado con nodos de Chebyshev en dirección z. 37 Ilustración 25 – Matrices de derivadas para nodos basados en DF y en Chebyshev. 38 Ilustración 26 – Mallado del dominio completo en el código sencillo. 38 Ilustración 27 – Mallado en Z según una función matemática. 39 Ilustración 28 – Mallado del dominio completo en el código complejo. 39 Ilustración 29 – Ejemplo de numeración de nodos con un índice [18]. 43 Ilustración 30 – Condiciones de contorno para el cálculo del potencial. 45 Ilustración 31 – Nodos que contienen a los electrodos. 50 Ilustración 32 – Condiciones de contorno adimensionales impuestas sobre los electrodos 50 Ilustración 33 – Nodos que contienen al electrodo visto. 54 Ilustración 34 – Condiciones de contorno adimensionales impuestas sobre el electrodo visto. 54 Ilustración 35 – Condiciones de contorno para resolver el problema fluido. 55 Ilustración 36 – Contorno y subcontorno del dominio. 57 Ilustración 37 – Densidad de carga obtenida con el código de Matlab. Código sencillo. 78 Ilustración 38 - Densidad de carga obtenida con el código de Matlab. Código complejo. 78 Ilustración 39 – Densidad de carga obtenida en [15] usando GHOST (Código desarrollado por la Universidad de Kentucky). 78 Ilustración 40 – Potencial eléctrico obtenido usando Matlab. Código sencillo. 79 Ilustración 41 – Potencial eléctrico obtenido usando Matlab. Código complejo. 80 Ilustración 42 – Potencial eléctrico obtenido en Suzen et Al [15] con GHOST. 80 Ilustración 43 – Potencial eléctrico obtenido en [20] con Oracle 81 Ilustración 44 – Potencial eléctrico en todo el dominio. Código sencillo. 81 Ilustración 45 – Potencial eléctrico en todo el dominio. Código complejo. 82 Ilustración 46 – Representación vectorial de la fuerza electrohidrodinámica entre los electrodos. Obtenida con el código de Matlab. 82 Ilustración 47 – Representación vectorial de la fuerza electrohidrodinámica. Publicada en [24]. 83 Ilustración 48 - Relación entre la fuerza y la distancia normal a la pared al comienzo del electrodo embebido. 83 Ilustración 49 – Relación entre la fuerza y la distancia normal a la pared al comienzo del electrodo embebido, publicado en [24]. 84 Ilustración 50 – Módulo de la fuerza en el entorno del actuador. 84 Ilustración 51 – Componentes x y z de la fuerza electrohidrodinámica. 85 Ilustración 52 – Región para el cálculo de la fuerza media. 86 Ilustración 53 – Fuerza máxima en función del espesor de dieléctrico. 86 Ilustración 54 – Fuerza máxima en z en función del espesor de dieléctrico. 87 Ilustración 55 – Fuerza media en x en función del espesor del dieléctrico. 87 Ilustración 56 – Fuerza máxima en función del gap. 88 Ilustración 57 – Fuerza media en x en función del gap. 88 Ilustración 58 – Relación fuerza máxima – voltaje. 89 Ilustración 59 – Relación velocidad inducida – voltaje [9]. 89 Ilustración 60 – Fenómeno de apantallamiento. 90 Ilustración 61 – Densidad de carga adimensional en función de la longitud de Debye. 91 Ilustración 62 – Densidad de carga máxima adimensional en función de la longitud de Debye. 92 Ilustración 63 – Fuerza adimensional en función de la longitud de Debye. 94 Ilustración 64 – Fuerza adimensional máxima en función de la longitud de Debye. 94 Ilustración 65 – Relación velocidad inducida – longitud de Debye [9]. 96 Ilustración 66 – Relación velocidad inducida – densidad de carga máxima, [9]. 96 Ilustración 67 – Problema de densidad de carga con condiciones de contorno periódicas. 97 Ilustración 68 – Problema potencial con condiciones de contorno periódicas. 97 18 Notación 𝑛𝑒 Densidad de electrones o densidad de carga 𝑛𝑖 Densidad de iones positivos 𝑛0 Densidad de plasma 𝑇𝑒 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 𝜈𝑒,𝑖 Frecuencia de colisiones entre iones y electrones 𝜈𝑒,𝑛 Frecuencia de colisiones entre neutrones y electrones 𝐷𝐵𝐷 Dielectric Barrier Discharge 𝑆𝐷𝐵𝐷 Single Dielectric Barrier Discharge 𝜌𝑐 Densidad de carga volumétrica 𝜌𝑐,max Densidad de carga volumétrica máxima 𝐸 󰇍 Campo eléctrico Φ Campo potencial 𝜀 Permitividad 𝜀𝑟 Permitividad relativa del medio 𝜙 Potencial debido a campos eléctricos externos 𝜑 Potencial debido a la densidad de carga en el propio 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 Ψ Función de corriente 𝑃𝑟 Número de Prandtl 𝐵 󰇍 Campo magnético 𝐷 󰇍 󰇍 Campo eléctrico inducido 𝐻 󰇍 󰇍 Campo magnético inducido 𝐽𝑓 󰇍 󰇍 󰇍 Corriente inducida 𝑤 Frecuencia de la corriente AC 𝜏 Tiempo adimensional 𝑢 Velocidad adimensional en dirección x 𝑣 Velocidad adimensional en dirección z 𝜌𝑠 Densidad de carga superficial ESTRUCTURA DEL TRABAJO El estudio está dividido en cinco secciones, el contenido de cada una de los bloques queda recogido a continuación:  Sección 1: Fundamentos teóricos. Se busca orientar al lector acerca de los conceptos básicos del plasma y del problema de control de flujo con actuadores de plasma. Se divide en subsecciones: 1.1. Conceptos básicos del plasma, 1.2. Control de flujo con generación de plasma. Física del proceso, 1.3. Principios, arquitectura, ventajas y limitaciones del actuador de plasma.  Sección 2: Modelado del problema. Se presentan las ecuaciones que gobiernan el problema, tanto fluido como eléctrico, así como la metodología empleada a lo largo del desarrollo del código numérico para simula la física del problema y sobre el que se trabaja para la obtención de resultados. Seguidamente, se trata el mallado realizado para discretizar el dominio, destacando sus principales características. Se describen los recursos matemáticos utilizados para la resolución numérica y se particularizan para cada problema. Por otro lado, se presentan las condiciones de contorno periódicas, que permitirán en el estudio del actuador con un caudal fluido, además de junto a infinitos actuadores adyacentes. Finalmente, se presentan los códigos elaborados junto con una breve explicación de su funcionamiento. La sección se divide en subsecciones: 2.1. Formulación del problema completo, 2.2. Mallado del dominio, 2.3. Resolución numérica del problema, 2.4. Particularización al problema eléctrico, 2.5. Particularización al problema fluido, 2.6. Problema con condiciones de contorno periódicas y 2.7. Descripción de los códigos elaborados.  Sección 3: Análisis de resultados. Se validan los resultados obtenidos al compararlo con modelos numéricos presentes en la literatura. Posteriormente se analiza la influencia de distintos parámetros de diseño de los actuadores DBD sobre su actuación. Se presentan además los resultados arrojados por el código al imponer condiciones de contorno periódicas. Por último, se estima el coste computacional del método elaborado. Se divide en las siguientes subsecciones: 3.1. Validación del código, 3.2. Influencia de distintos parámetros sobre el comportamiento del actuador, 3.3. Influencia de los parámetros experimentales del modelo SH sobre el comportamiento del actuador, 3.4. Problema eléctrico con condiciones de contorno periódicas, 3.5. Coste computacional.  Sección 4: Conclusiones. En esta sección se reflejan las conclusiones a las que se han llegado a lo largo de la investigación y se comparan con los objetivos iniciales propuestos al comienzo del estudio.  Sección 5: Líneas futuras. Finalmente, se tratan las posibilidades que quedan abiertas al estudio tras la investigación y se detallan ciertas vías de ampliación del trabajo llevado a cabo. 12 12 ESTADO DEL ARTE Nothing in life is to be feared, it is only to be understood. Now is the time to understand more, so that we may fear less. Marie Curie, 1973. En las últimas décadas, la investigación y desarrollo de dispositivos de control de flujo poco intrusivos ha experimentado un notable crecimiento, tratando e involucrando ámbitos como la acústica o la generación de plasma, en la cual se basa este estudio. El control de la capa límite con actuadores de plasma es uno de los temas de investigación del momento, desde que en 1998 empezaran a publicarse los primeros escritos [1,2]. Debido a la simplicidad del dispositivo, así como el potencial que se augura, este tipo de actuadores ha ido adquiriendo cada vez más interés en el ámbito de la ingeniería aeronáutica. El actuador de plasma DBD (Dieléctric Barrier Discharge), es capaz de actuar sobre el flujo, induciendo componentes de velocidad sobre este. Entre las múltiples aplicaciones recogidas en la literatura es capaz de retrasar el punto de desprendimiento de la capa límite en un perfil [3,4], lo cual mejora la sustentación global de este, reducir la resistencia aerodinámica en un cuerpo romo [5,6], disminuir los niveles de ruido [7] o eliminar fenómenos como las ondas de Tollmien-Schlichting [8]. Además, se caracterizan por ser dispositivos ligeros, compactos y de arquitectura extremadamente sencilla, lo que los hace robustos. En la actualidad, se está investigando la implementación de este tipo de actuadores como sustitutos de los sistemas de control de vuelo tradicionales, es decir, superficies de control móviles de las aeronaves, con la considerable reducción del peso que eso conlleva, así como el aumento del espacio para la acumulación de combustible, lo que supondría el incremento del alcance del vehículo. También se ha depositado gran esperanza en el control de la capa límite de álabes de turbina, lo que permitiría aumentar la eficiencia de los motores de reacción, así como la implementación de estrategias de control activo que expandan los límites de actuación de la planta de potencia. La lista de aplicaciones de estos dispositivos es interminable, siendo capaces de retrasar la ruptura de la barrera del sonido en un avión supersónico o incluso de ayudar a la disipación de calor sobre ciertos puntos del fuselaje de las aeronaves. Sin embargo, los modelos numéricos para el estudio de la física del fenómeno de generación de plasma y la influencia de este sobre un fluido pueden ser extraordinariamente complejos y requerir un gran esfuerzo computacional, como es el caso de los códigos con modelos de plasma microscópicos, que emulan con detalle la dinámica del proceso de las distintas especies involucradas. En 2005, se introdujo el modelo de Suzen- Huang (SH), un método numérico macroscópico, mucho más sencillo y que permite obtener resultados interesantes desde el punto de vista ingenieril, con muchos menos recursos. El código para la simulación del actuador elaborado en este trabajo se basa en dicho modelo. 14 14 MOTIVACIÓN DEL PROBLEMA El control de flujo mediante dispositivos mecánicos, como podría ser la deflexión de una corriente fluida con un perfil aerodinámico o el soplado en la capa límite, es un problema ampliamente conocido; una metodología que presenta inconvenientes como la complejidad de los mecanismos, su elevado coste o el tamaño de los dispositivos de control, que pueden resultar intrusivos para la corriente y reducir la eficiencia del proceso, así como su limitado rango de actuación. Estos son los aspectos más desfavorables y cuya mejora resultaría de gran interés para innumerables aplicaciones aeronáuticas. La utilización del plasma implica la completa desaparición de los dispositivos de control mecánico, actuando sobre la corriente a través de fuerzas electromagnéticas. Esta tecnología se implementa por medio de actuadores de arquitectura muy sencilla, sin partes móviles y que no obstaculizan al movimiento del fluido, lo que se traduce en reducción de la resistencia y el consiguiente incremento de la eficiencia en el proceso de control. Recientemente Corke, un reconocido investigador en el ámbito de los actuadores de plasma, realizó un estudio numérico donde implementaba este tipo de dispositivos sobre un UAV Predator con el fin de reducir la resistencia aerodinámica, el resultado fue la disminución del 65% de la resistencia, así como el incremento del alcance de la aeronave en un 300% y de la autonomía, que paso de 24 a 36 horas. Resulta indudable que el actuador de plasma presenta un gran potencial en infinidad de ámbitos y que está en ferviente desarrollo. Es por ello que existen complejos códigos y softwares que permiten simular el comportamiento de estos dispositivos, pero la mayoría de ellos requieren gran capacidad computacional y largos tiempos de cálculo. El objetivo de este trabajo es desarrollar un método de simulación que permita, en primera aproximación, obtener resultados ingenieriles que den idea del comportamiento del actuador ante distintas condiciones, y de una manera sencilla, sin necesidad de grandes recursos computacionales, permitiendo obtener resultados aproximados de forma rápida sin más que un ordenador personal y el software Matlab. Para ello se elabora un modelo numérico simple que, a diferencia de los códigos tipo “caja negra” convencionales, facilite la investigación a aquellos que no dispongan de la capacidad computacional que los métodos comerciales requieren. Además, se busca la creación de un modelo flexible y fácil de comprender, sobre el que el usuario pueda modificar los parámetros, ecuaciones y condiciones de contorno a su antojo. Debido a la extensión del proyecto el modelado del actuador DBD se ha dividido en dos trabajos de fin de grado. En este proyecto se centra especialmente en el estudio, implementación y validación del problema eléctrico, necesario para obtener la solicitación generada por el actuador sobre el aire en el entorno de este. La segunda parte del modelo ha sido elaborada por Pablo Dugo Domínguez, uno de mis mejores amigos, encargado de la implementación del problema fluido y del estudio del campo de velocidades inducidas por el actuador, tarea que queda recogida en su trabajo de fin de grado [9]. 22 22 corriente.  Arco Eléctrico: Esta última zona corresponde a descargas de alta corriente, superior a 100 A, y que son muy brillantes. En el diagrama anterior hay dos valores de tensión característicos:  Tensión de ignición de gases: Es la tensión a partir de la cual la descarga se vuelve visible, es decir, la tensión en el punto D del diagrama anterior.  Tensión de ruptura de gases: Es la tensión a partir de la cual se producen arcos eléctricos, es decir, la tensión en el punto H del diagrama voltaje-intensidad. Plasmas de Descarga de Barrera Dieléctrica La descarga de Barrera Dieléctrica es una descarga no térmica producida al aplicar alto voltaje entre dos electrodos separados entre sí por al menos una pequeña capa de material dieléctrico, que evita que la descarga se transforme en un arco eléctrico. Debido a las capas de material aislante la descarga se realiza con voltajes alternos, provocando una descarga cada semi ciclo del voltaje, ya que una corriente continua no puede atravesar el material dieléctrico. Las capas aislantes actúan como limitadores de la corriente, evitando que se formen arcos eléctricos y permitiendo la formación de plasma en el espacio libre entre los electrodos (gap). Este tipo de descargas son posibles en un amplio rango de presiones, desde ~ 0,1 bar hasta ~ 10 bar, y temperaturas. Resultando de especial interés, debido a sus numerosas aplicaciones industriales, las descargas a presión atmosféricas. Las descargas DBD se pueden clasificar en dos modos de operación: la descarga filamentosa y la difusa luminiscente. El primer modo se obtiene a presión atmosférica, el plasma está formado por micro-descargas con canales conductores de ~ 0,1 mm y tiene duración de nanosegundos. Ilustración 6 – Descarga de Barrera Dieléctrica Filamentosa [28]. Por el contrario, el modo luminiscente se produce a baja presión y todo el espacio entre los electrodos está lleno de plasma. Ilustración 7 – Descarga de Barrera Dieléctrica Difusa. Los actuadores de plasma pueden tener distintas configuraciones, dependiendo de si trabajan volumétrica o superficialmente. En este estudio se tratará el segundo caso, actuadores superficiales, que generen la descarga sobre la propia superficie sobre la que se desea controlar el fluido. Ilustración 8 – Actuadores DBD superficiales, con descarga única y doble. Materiales dieléctricos Un dieléctrico es un aislante, es decir un material con baja conductividad eléctrica, que se puede polarizar al aplicar un campo eléctrico sobre él. Es importante mencionar que todos los materiales dieléctricos son aislantes, sin embargo, no todos los aislantes son dieléctricos. La diferencia entre ambos reside en la facilidad para polarizar a los primeros, propiedad que se puede cuantificar mediante la permitividad relativa, llamada también constante dieléctrica del material. Cuando el dieléctrico está bajo la acción de un campo eléctrico las cargas no viajan a través del material como en el caso de los conductores, sino que pasan de la posición de equilibrio medio a una posición polarizada, donde las cargas positivas se orientan en la dirección del campo eléctrico y las negativas en la dirección opuesta. Esto genera un campo eléctrico interno de sentido contrario, y el campo eléctrico total resultante es de menor intensidad. 24 24 Ilustración 9 – Material dieléctrico sin y con presencia de un campo eléctrico. En este estudio se utilizará el Kapton como material dieléctrico entre los electrodos, sin bien podría sustituirse en el modelo sin más que cambiar su constante característica. En la siguiente tabla aparecen una serie de materiales y sus correspondientes permitividades relativas. Material Permitividad relativa Mica 6.30-9.30 Vidrio 4.90-7.50 Baquelita 3.50-6.00 Madera (Roble) 3.30 Kapton 2.70-3.50 Teflon 2.15 Tabla 1 – Constantes dieléctricas de distintos materiales. Fundamentos del actuador El control del flujo se basa en modificar el comportamiento del fluido en busca de una nueva configuración, para ello existen dispositivos pasivos de control de corriente y dispositivos activos. Los primeros son dispositivos típicamente mecánicos y que no añaden energía al sistema, mientras que los segundos sí que añaden energía al sistema por medio de un actuador, que a su vez necesita alimentación. Los actuadores de plasma tipo DBD son dispositivos de control activo, que como su propio nombre indica, se basan en la generación de plasmas no térmicos a presión atmosférica. Para su funcionamiento se utiliza corriente alterna, las descargas aparecen en la superficie del dieléctrico sobre el electrodo embebido, produciendo una inyección de momento, proveniente de las colisiones entre las partículas cargadas del plasma y las partículas neutras del aire, en el entorno del actuador en cada semi-ciclo de la corriente aplicada. Aunque la inyección de momento ocurre de forma alterna, el momento medio se dirige desde el electrodo expuesto al electrodo enterrado, derivando en la eyección del aire aguas abajo, alcanzando velocidades en torno a ~ 1-10 m/s y un momento medio del orden de ~ 0.01 – 0.1 N/m. En numerosos estudios se ha demostrado la viabilidad de estos dispositivos como elementos de control de flujo. Ejemplo de ello son las investigaciones llevadas a cabo por Corke et al. [10] donde se analiza el aumento de sustentación en una sección de ala, o el control de la separación de la capa límite en los álabes de turbinas de baja presión, estudiado por Hultgren y Ashpis [11]. Los resultados obtenidos en los diversos experimentos llevados a cabo desde la aparición de este tipo de actuadores demuestran que existen ciertos parámetros fundamentales para el control efectivo del flujo. Dichos parámetros son: - Posición, tamaño y orientación de los electrodos. - Hueco (gap) entre los electrodos. - Voltaje y frecuencia que aplica la fuente de corriente alterna. Por tanto, el control de flujo deriva en un problema complejo, con un gran número de variables de control que optimizar. Arquitectura del actuador Los actuadores con generación de plasma resultan ser dispositivos extremadamente sencillos, compuestos fundamentalmente por dos electrodos separados entre sí por un dieléctrico y sometidos a una diferencia de potencial que varía en el tiempo. La alimentación de este tipo de actuador es mediante corriente alterna, con frecuencias de entre ~ 1-100 kHz y voltajes del orden de 10 kV. Existen infinidad de configuraciones posibles, con electrodos de distintos tamaño y orientaciones, incluso se implementan arrays de actuadores simples, formando dispositivos de control más sofisticados. En lo que sigue, se elige el dispositivo de control más sencillo, denominado actuador SDBD (‘Single Dielectric Barrier Discharge’). Se compone de un electrodo visto, en contacto con el aire, un dieléctrico y un electrodo embebido entre el dieléctrico y la superficie sobre la que se desea controlar el flujo. Ilustración 10 – Arquitectura básica del actuador. Ventajas del actuador DBD Los actuadores aerodinámicos en general, ya bien pasivos o activos, buscan aumentar la sustentación, reducir la resistencia, retrasar el desprendimiento o disminuir los niveles de ruido entre otras funciones. Los actuadores DBD tienen un gran potencial, consiguiendo realizar las funciones anteriormente expuestas, pero de una forma más eficiente, con menor impacto medioambiental y mayor ahorro de energía y combustible. Para ello, se basan en la inyección de momento en el aire cercano a la superficie sobre la que se sitúa el dispositivo, siendo capaz de modificar el estado del flujo y controlar la corriente. Las principales ventajas respecto a los dispositivos tradicionales de control de flujo residen en su tamaño reducido, estando prácticamente incrustado sobre la superficie aerodinámica de forma que no perturba a la corriente, en su bajo peso y en la ausencia de partes móviles, que lo convierte en un dispositivo robusto, fiable y barato. Además, utilizan estrategias de control gobernadas por un controlador, lo que les permite dar una respuesta rápida, adaptarse a cambios fácilmente y la posibilidad de implementar aerodinámica activa, con leyes de control a tiempo real. Los actuadores DBD son dispositivos de especial interés para las generaciones venideras en la industria aeronáutica, posicionándose como un fuerte candidato para las funciones de control de la corriente en los futuros More Electric Aircrafts. 26 26 Limitaciones del actuador DBD Los actuadores DBD siguen en desarrollo, y, por tanto, siguen presentando limitaciones que impiden su implementación a gran escala en la industria. Una de las principales desventajas es la necesidad de pesados generadores de corriente con alto voltaje y la baja eficiencia que presentan los actuadores para convertir la energía eléctrica en energía cinética. Sin embargo, la evolución hacia los More Electric Aircraft supondrá la renovación de los medios de generación de corriente en las aeronaves, beneficiando a la implementacion de este tipo de dispositivos. Otra limitación es la incapacidad de producir con un actuador simple chorros de aire inducido suficientemente rápidos como para satisfacer las necesidades de las aplicaciones aeronáuticas (~ 80 m/s), lo que obliga a utilizar matrices con múltiples actuadores. La mayor parte del esfuerzo de investigación en la actualidad reside en conseguir dispositivos que permitan convertir la energía de forma más eficiente y la optimización de los arrays de actuadores, apareciendo nuevas soluciones de control como los canales de flujo impulsado con plasma, que pueden ser usados para crear succión o eyectar corriente en aplicaciones de control de flujo. Ilustración 11 – Canal tipo DBD de flujo impulsado por plasma. 2 MODELADO DEL PROBLEMA El modelado de un actuador DBD se puede abordar desde dos perspectivas diferentes, la microscópica, donde destacan las publicaciones de Massines et al. (1998) [12], Boeuf et al. (2005) [13] y Rouffet (2008) [14], y la macroscópica, donde resultan relevantes los estudios llevados a cabo por Roth et al. (1998) [1,2] y Suzen et al (2005) [15]. En el primer caso se centra en conocer el movimiento de las partículas cargadas, calculando la trayectoria seguida por cada una ellas. Se centra en procesos definidos por una escala temporal del orden de 10−8−10−9 s, que concuerda con la escala de tiempo de la movilidad de las partículas, y los cálculos requieren una gran capacidad computacional. Por el contrario, los modelos macroscópicos no se basan en un estudio tan detallado de la física del problema, sino que apuntan a datos sobre el flujo de corriente inducido por el actuador. En el caso de Suzen et Al. [15] para estudiar el efecto de la densidad de carga sobre el electrodo embebido se considera esta como un único ente y no como un conjunto de partículas, cada una con un comportamiendo determinado. La escala temporal de este tipo de métodos es del orden de 10−4 s, quedando definida por la frecuencia de la corriente AC aplicada en los electrodos, y requieren menos recursos para la resolución del problema. En este proyecto se opta por elaborar un método macroscópico sencillo, para el cual será necesario el cálculo de las magnitudes asociadas al problema eléctrico y de las ecuaciones de Navier-Stokes, que gobiernan el comportamiento del fluido. En la primera parte de esta sección se trata la formulación y modelado del problema. Posteriormente, se presenta el mallado utilizado para discretizar el dominio, y finalmente se detallan las herramientas matemáticas necesarias para la resolución numéricar de las ecuaciones, tanto del problema eléctrico como del fluido. Por último, se incluyen los códigos elaborados con una breve explicación de su funcionamiento. 2.1 Formulación del problema completo El estudio del comportamiento del plasma puede ser un problema extremadamente complejo, con infinidad de comportamiento distintos que hacen que su modelado se convierta en una tarea ardua y sólo abarcable por medio de sistemas complejos, con modelos multi-especie. Estos sistemas lindan entre los comportamientos ordenados y los caóticos, si bien la física del plasma en la mayoría de los casos no puede ser descrita mediante las ecuaciones fluidas simples, tampoco se rige por una absoluta aleatoriedad. Sin embargo, el objetivo de este trabajo no es un estudio pormenorizado de la dinámica del cuarto estado de la materia, la investigación se centra en el modelado, de forma sencilla, de un actuador DBD. Por tanto, en lo que sigue se desarrolla un modelo de actuador de plasma simple inspirado en el modelo de Suzen-Huang [15] (SH en adelante). El modelo SH fue publicado en 2005 y representa el comportamiento de un actuador DBD, permitiendo obtener resultados interesantes desde el punto de vista ingenieril con pocos recursos y esfuerzo computacional. Está basado en el flujo electroosmótico, un concepto introducido por Henry en 1931 [16], y que se asocia al movimiento inducido sobre un fluido encerrado en un tubo capilar sobre el que se aplica cierto voltaje. En los actuadores de plasma, donde la longitud de Debye es pequeña, el movimiento fluido se asemeja al del flujo electroosmótico. La principal característica del modelo SH es que descompone el potencial del que deriva el campo eléctrico en dos términos, el asociado al campo eléctrico externo que se aplica sobre los electrodos y el potencial debido a las cargas que se acumulan sobre la superficie del actuador. Así, se obtienen las dos ecuaciones características del problema eléctrico, que permiten posteriormente calcular la fuerza volumétrica generada por el actuador y que se incluye finalmente en las ecuaciones de Navier-Stokes como término fuente, permitiendo conocer el comportamiento del fluido en presencia de dicha solicitación. A continuación, se lleva a cabo el desarrollo y justificación de las ecuaciones que gobiernan el comportamiento del actuador, así como de las condiciones de contorno impuestas. Para modelar el efecto del actuador sobre el fluido se deben resolver las ecuaciones incompresibles de Navier- Stokes, las cuales aparecen a continuación: 28 28 ∇·𝑣 =0 (2) 𝜌𝜕𝑣 𝜕𝑡+𝜌(𝑣 ·∇𝑣 )−𝜇∇2𝑣 +∇𝑝=𝐹𝑏 󰇍 󰇍 󰇍 󰇍 (3) Donde 𝜌 es la densidad del aire [𝑘𝑔 𝑚3], 𝑣 la velocidad [𝑚/𝑠], 𝜇 la viscosidad dinámica del aire [𝑘𝑔 𝑚·𝑠], 𝑝 la presión [𝑃𝑎] y 𝐹𝑏 󰇍 󰇍 󰇍 󰇍 es la fuerza volumétrica de Lorentz [𝑁 𝑚3], la cual se puede obtener mediante la ecuación de Lorentz. 𝐹𝑏 󰇍 󰇍 󰇍 󰇍 =𝜌𝑐(𝐸 󰇍 +𝑣 ×𝐵 󰇍 ) (4) Siendo 𝜌𝑐 la densidad de carga [𝐶 𝑚3], 𝐸 󰇍 el campo eléctrico [𝑉 𝑚] y 𝐵 󰇍 el campo magnético [𝑇]. La expresión de la fuerza volumétrica se puede desarrollar aun más recurriendo a las ecuaciones de Maxwell en forma diferencial: ∇·𝐷 󰇍 󰇍 =𝜌𝑐 (5) ∇·𝐵 󰇍 =0 (6) ∇×𝐸 󰇍 =−𝜕𝐵 󰇍 𝜕𝑡 (7) ∇×𝐻 󰇍 󰇍 =𝐽 𝑓+𝜕𝐷 󰇍 󰇍 𝜕𝑡 (8) Donde 𝐷 󰇍 󰇍 es el campo eléctrico inducido [𝑉 𝑚], 𝐻 󰇍 󰇍 es el campo magnético inducido [𝑇] y 𝐽 𝑓 la corriente inducida [𝐴 𝑚]. En el caso de la descarga de gases el sistema se considera cuasi-estacionario, ya que el periodo de la corriente AC es mucho mayor (10−4 s) que el tiempo que necesitan las cargas para redistribuirse en el plasma (10−8−10−9 s). Esta hipótesis implica que 𝐽 𝑓,𝐵 󰇍 y 𝐻 󰇍 󰇍 son nulos, y además, que las derivadas temporales también lo son. Por tanto, sólo se mantienen las ecuaciónes (5) y (7). El campo eléctrico y el campo eléctrico inducido están relacionador por medio de la expresión: 𝐷 󰇍 󰇍 =𝜀𝐸 󰇍 (9) Siendo 𝜀 el coeficiente dieléctrico [𝐹 𝑚], que se puede expresar como el producto de la permitividad dieléctrica relativa del medio por la del vacío, 𝜀=𝜀𝑟·𝜀0. Además, la ecuación (7) quedaría de la siguiente forma: ∇×𝐸 󰇍 =0 (10) Por tanto, se puede concluir que el campo eléctrico deriva de un potencial: 𝐸 󰇍 =−∇Φ (11) Lo cual permite simplificar la ecuación de la fuerza de Lorentz: 𝐹𝑏 󰇍 󰇍 󰇍 󰇍 =𝜌𝑐𝐸 󰇍 =−𝜌𝑐∇Φ (12) Donde Φ representa al potencial eléctrico. Si el plasma generado tiene una grado de ionización bajo (𝛼≪1), es decir, plasma débilmente ionizado; se puede descomponer el potencial en dos partes: el debido a campos eléctricos externos, 𝜙, y potencial debido a la densidad de cargas sobre la superficie del dieléctrico, 𝜑. Φ=𝜙+𝜑 (13) El modelo SH tiene en cuenta la aproximación de Debye-Huckel para el tratamiento de la doble capa eléctrica, que es la región de interfase entre el medio fluido y el dieléctrico. Esta aproximación sugiere que, para longitudes de Debye pequeñas, el potencial debido a las cargas concentradas en la superficie del dieléctrico, 𝜑, sólo influye en la distribución de las propias partículas cargadas y no en el campo eléctrico sobre el dominio. Como en el caso del actuador de plasma, la longitud de Debye y la carga sobre la superficie del dieléctrico son pequeñas, el potencial eléctrico total se divide en dos términos que actúan de forma independiente sobre la densidad de carga y sobre el campo eléctrico. Por tanto, introduciendo la ecuación (13) en (5) se obtienen las principales ecuaciones del modelo SH. ∇·(𝜀𝑟∇𝜙)=0 (14) ∇·(𝜀𝑟∇𝜑)=−𝜌𝑐 𝜀0 (15) Se puede comprobar, tal y como se mencionó con anterioridad, que el potencial 𝜑 está asociado a la densidad de carga sobre la superficie del dieléctrico, 𝜌𝑐, y que el campo eléctrico pasa a depender exclusivamente de 𝜙, lo que a su vez permite simplificar la expresión de la fuerza de Lorentz. 𝐸 󰇍 =−∇𝜙 (16) 𝐹𝑏 󰇍 󰇍 󰇍 󰇍 =𝜌𝑐·(−∇𝜙) (17) Resolviendo la ecuación (14) se obtiene el potencial debido al campo eléctrico aplicado sobre los electrodos, el cual es necesario para el posterior cálculo de la fuerza electrohidrodinámica, para su resolución se aplican las siguientes condiciones de contorno: Ilustración 12 – Dominio y condiciones de contorno para resolver la ecuación (14). Cabe destacar que durante el trabajo se realizarán dos versiones del código, una más sencilla donde no se tendrá en cuenta el efecto del espesor de los electrodos sobre el problema y otra que sí tendrá en cuenta el espesor de los electrodos, así como un mallado más fino. Por otro lado, se puede reescribir la ecuación (15) exclusivamente en función de la densidad de carga. Si consideramos el caso en el que durante la descarga sólo estén presentes, de forma significativa, protones y Aire: ∇ (𝜀𝑟1·∇𝜙 )=0 𝜀𝑟1=1 Contorno exterior: ∂ϕ ∂n=0 𝜙=𝜙(𝑡) 𝜙=0 𝜀𝑟=𝑚𝑒𝑑𝑖𝑎(𝜀𝑟1,𝜀𝑟2) Dieléctrico (Kapton): ∇ (𝜀𝑟2·∇𝜙 )=0 𝜀𝑟2=2.7 30 30 electrones podemos aproximar la densidad de carga neta, 𝜌𝑐, dentro del plasma según la ley estadística de Boltzmann para el equilibrio térmico entre electrones y protones. 𝑛𝑒=𝑛0·exp(𝑒𝜑 𝑘𝑇𝑒) (18) 𝑛𝑖=𝑛0·exp(−𝑒𝜑 𝑘𝑇𝑖) (19) 𝜌𝑐=𝑒(𝑛𝑖−𝑛𝑒) (20) 𝜌𝑐=𝑒𝑛0(exp(−𝑒𝜑 𝑘𝑇𝑖)−exp(𝑒𝜑 𝑘𝑇𝑒)) (21) Donde 𝑛𝑖 es la densidad de iones positivos, 𝑛𝑒, la densidad de electrones y 𝑛0 la densidad de plasma. Aplicando el desarrollo en serie de Taylor para la función exponencial se llega a: 𝜌𝑐=−𝑒2𝑛0[1 𝑘𝑇𝑖+1 𝑘𝑇𝑒] 𝜑 (22) Se introduce ahora la expresión de la longitud de Debye (𝜆𝐷), que como se definió anteriormente, es la distancia que determina el radio de la esfera de influencia de cada carga dentro del plasma: 1 𝜆𝐷=𝑒2𝑛0 𝜀0[1 𝑘𝑇𝑖+1 𝑘𝑇𝑒] (23) Combinando ambas expresiones, obtenemos la expresión del potencial eléctrico debido a la densidad de carga en el plasma en función de la longitud de Debye. 𝜑=−𝜌𝐶𝜆𝐷 2 𝜀0 (24) Introduciendo la expresión (24) en la ecuación (15) se llega a la ecuación del modelo SH para la distribución de densidad de carga: ∇·(𝜀𝑟∇𝜌𝑐)=𝜌𝑐 𝜆𝐷 2 (25) Esta expresión permite conocer la densidad de carga sobre la superficie del dieléctrico. Las condiciones de contorno impuestas en este caso vienen dadas por: Ilustración 13 – Dominio y condiciones de contorno para calcular la densidad de carga. Las condiciones de contorno impuestas son ligeramente distintas a las presentadas por Suzen et Al. en [15], se ha modificado la condición establecida sobre el electrodo visto, considerando que la densidad de carga en el interior de este es nula. Sobre la superficie de dieléctrico que cubre al electrodo embebido el modelo establece una distribución de densidad de carga basada en la experimentación, que varía a lo largo del electrodo y con el tiempo. 𝜌𝑐,𝑤(𝑥,𝑡)=𝜌𝑐𝑚𝑎𝑥 𝐺(𝑥) 𝑓(𝑡) (26) Donde 𝜌𝑐𝑚𝑎𝑥 es el máximo valor de densidad de carga sobre el dominio, un parámetro determinado experimentalmente, G(x) es una distribución semi Gaussiana que modela la distribución del plasma a lo largo del electrodo embebido y 𝑓(𝑡) es la variación temporal del voltaje, que dependerá del tipo de onda que se aplique, en este caso se consideran ondas sinusoidales pero bien podrían ser ondas rectangulares u otros tipos, su expresión sería: 𝑓(𝑡)=sin (2𝜋𝑤𝑡) (27) Siendo 𝑤 la frecuencia de la corriente AC aplicada [Hz]. En cuanto a la distribución semi Gaussiana, se obtiene al analizar la emisión de luz del plasma al producirse la descarga, esta es proporcional a la densidad de carga y los resultados obtenidos se asemejan con la distribución 𝐺(𝑥) impuesta en el modelo SH. Ilustración 14 – Medida 3D de la emisión de luz durante un ciclo completo de corriente AC [17]. Aire: ∇ (𝜀𝑟1 ∇𝜌𝑐 )=𝜌𝑐 𝜆𝑑 2 𝜀𝑟1=1 Contorno exterior: 𝜌𝑐=0 𝜕𝜌𝑐 𝜕𝑛=0 𝜕𝜌𝑐 𝜕𝑛=0 𝜌𝑐=0 𝜌𝑐=𝜌𝑐,𝑤(𝑥,𝑡) Dieléctrico (Kapton): ∇ (𝜀𝑟2·∇𝜌𝑐)=0 𝜀𝑟2=2.7 38 38 Los nodos son los ceros del polinomio y, por tanto: 𝑧𝑗=𝑎+𝑏 2+𝑎−𝑏 2·cos[𝜋 (𝑗−1) 𝑁−1] 𝑑𝑜𝑛𝑑𝑒 𝑗=1…𝑁𝑧 𝑦 𝑧𝑗∈[𝑎,𝑏] (47) El principal inconveniente del mallado con Chebyshev es que las matrices de derivadas son llenas, lo que limita el número de puntos que se pueden establecer a lo largo del eje Z y con ello la precisión del código, puesto que un ligero aumento en el número de nodos implica el llenado de la memoria del ordenador y que no pueda ejecutarse el programa. Es por ello, que en la versión más compleja se haya optado por la utilización de diferencias finitas a lo largo del eje Z, puesto que en este caso las matrices de derivadas son vacias y permite aunar una gran cantidad de nodos, incrementando la precisión y fiabilidad de los resultados a costa de aumentar también el esfuerzo computacional. Ilustración 25 – Matrices de derivadas para nodos basados en DF y en Chebyshev. Finalmente, el mallado completo, tendría un aspecto similar al mostrado a continuación, donde se representa una discretización con un menor número de celdas para facilitar su visualización. Ilustración 26 – Mallado del dominio completo en el código sencillo. Para el código más elaborado el mallado a lo largo del eje Z es más preciso. En la región fluida concentra mayor número de celdas en el entorno de los electrodos y cerca de la pared del dieléctrico, mientras que en la región dieléctrica presenta un mallado equiespaciado. Para crear la malla del medio fluido se elabora una trama que varía según la función matemática 𝑦=𝑥2, lo cual permite concentrar un mayor número de celdas en la cercanía de la pared y hace que la malla vaya creciendo a medida que se separa de esta. El fundamento de hacer que el mallado respete el comportamiento de dicha norma matemática se puede entender con la siguiente ilustración. Ilustración 27 – Mallado en Z según una función matemática. Como se puede apreciar, más de la mitad de los nodos establecidos se concentran en las inmedianciones de la pared. Para implementar este mallado en Matlab se creado un vector equiespaciado entre el límite inferior y el superior de nuestro dominio y se ha proyectado sobre otro vector equiespaciado que en este caso va desde 0 a 1 y está elevado al cuadrado, siguiendo así la función matemática deseada. A continuación, se representa una imagen del mallado completo, al igual que antes reduciendo el número de celdas para facilitar la visualización del mismo. Ilustración 28 – Mallado del dominio completo en el código complejo. 40 40 2.3 Resolución numérica del problema Una vez establecidas las ecuaciones que rigen el comportamiento del actuador, así como la discretización del dominio sobre la que se llevarán a cabo los cálculos es necesario recoger las herramientas matemáticas que se emplean en la resolución del problema. En primer lugar, se presenta de manera general el método de colocación, que permite la resolución numérica de ecuaciones diferenciales ordinarias y ecuaciones en derivadas parciales entre otras. Posteriormente, se particulariza dicho método hasta adaptarlo a la formulación del problema eléctrico. Finalmente, se vuelve a adaptar el método de colocación junto con la función de corriente, esta vez para resolver el problema fluido. El método de colocación El método de colocación es un método numérico aplicado tanto a ecuaciones de evolución como ecuaciones estacionarias cuya solución es regular, es decir, sin “saltos bruscos”. Se basa en expresar las funciones incógnitas del problema como combinaciones lineales de funciones conocidas, como podrían ser los polinomios interpolantes de Lagrange. Por tanto, las ecuaciones incógnitas quedarían reescritas como expansiones de coeficientes desconocidos en un primer momento. Los coeficientes se determinan sustituyendo dichas expansiones y sus derivadas en las ecuaciones diferenciales y condiciones de contorno del problema e imponiendo que se verifiquen en un conjunto dado de puntos del dominio de integración. Dichos puntos, que son prefijados, se denominan puntos de colocación. El resultado es un sistema con tantas ecuaciones como incógnitas que permite calcular los coeficientes de la expansión y resolver el problema. Método de colocación para problemas de una variable con interpolantes de Lagrange y nodos de Chebyshev A pesar de que las ecuaciones que rigen el comportamiento del actuador de plasma tienen más de una variable, se introduce este problema como ayuda para la comprensión del funcionamiento del método. Considerese un problema diferencial de una variable independiente del que se desea conocer su solución en un intervalo [a, b], sean además f1, f2, …, fN los valores de la función incógnita en N puntos prefijados del intervalo, a = x1 < x2 < x3 < … < xN-1 < xN = b. Se busca obtener una aproximación, f(x), de la solución del problema en [a, b] mediante una expansión, 𝑓(𝑥)=∑𝐿𝑖(𝑥)·𝑓𝑖 𝑁 𝑖=1 (48) Donde 𝐿𝑖(𝑥) son polinomios interpolantes de Lagrange para cada uno de los puntos dados, 𝑖∈{1,…,𝑁}. Los polinomios de Lagrange son de la forma: 𝐿𝑖(𝑥)=∏𝑥−𝑥𝑗 𝑥𝑖−𝑥𝑗 𝑁 𝑗=1 𝑗≠𝑖 𝑖=1,…,𝑁 (49) Este tipo de polinomios se comporta según la delta de Kronecker, es decir, 𝐿𝑖(𝑥𝑗)=𝛿𝑖𝑗 ∀ 𝑖,𝑗=1,…,𝑁. 𝛿𝑖𝑗= {1 𝑠𝑖 𝑖=𝑗 0 𝑠𝑖 𝑖≠𝑗 (50) De la propia definición del polinomio de Lagrange se obtiene el valor de su derivada para cualquier nodo xi del intervalo, que viene dado a continuación. 𝐿𝑖′(𝑥)=∑(𝑥−𝑥1)…(𝑥−𝑥𝑖−1)(𝑥−𝑥𝑖+1)…(𝑥−𝑥𝑘−1)(𝑥−𝑥𝑘+1)…(𝑥−𝑥𝑁) (𝑥𝑖−𝑥1)…(𝑥𝑖−𝑥𝑖−1)(𝑥𝑖−𝑥𝑖+1)…(𝑥−𝑥𝑁) 𝑁 𝑘=1 (51) 𝑓𝑖′=∑𝑑𝐿𝑗(𝑥) 𝑑𝑥 | 𝑁 𝑗=1 𝑥=𝑥𝑖 𝑓𝑗′=∑𝐿𝑗′(𝑥𝑖) 𝑓𝑗 𝑖=1…𝑁 𝑁 𝑗=1 (52) Si se hace el paso análogo y se aproxima la derivada de la función incógina mediante un polinomio interpolador de Lagrange, g(x), se tiene que, 𝑔(𝑥)=∑𝐿𝑗(𝑥)·𝑓𝑗′ 𝑁 𝑗=1 (53) Que permite a su vez, aproximar los valores de la derivada segunda como: 𝑓𝑖′′=𝑔(𝑥𝑖)′=∑𝐿𝑗′(𝑥𝑖)·𝑓𝑗′= ∑𝐿𝑗′(𝑥𝑖)·(∑𝐿𝑘 ′(𝑥𝑗)·𝑓𝑗) 𝑁 𝑘=1 𝑁 𝑗=1 𝑁 𝑗=1 (54) El valor de las funciones y sus derivadas puede escribirse vectorialmente, lo que simplifica la notación y permite escribir las ecuaciones de forma compacta. 𝑓 =[𝑓1 𝑓2 … 𝑓𝑁], 𝑓 ′=[𝑓1′ 𝑓2′ … 𝑓3′], 𝑓 ′′=[𝑓1′′ 𝑓2′′ … 𝑓3′′] 𝑓 ′=𝐿𝑝𝑥·𝑓 𝑓 ′′=𝐿𝑝𝑥·𝑓 ′=𝐿𝑝𝑥 2·𝑓 𝑓 𝑛=𝐿𝑝𝑥·𝑓 𝑛−1=𝐿𝑝𝑥 𝑛·𝑓 (55) Siendo 𝐿𝑝𝑥 la matriz con las derivadas de los polinomios de Lagrange: 𝐿𝑝𝑥=(𝐿1′(𝑥1)⋯ 𝐿𝑁 ′(𝑥1) ⋮ ⋱ ⋮ 𝐿1′(𝑥𝑁)⋯ 𝐿𝑁 ′(𝑥𝑁)) (56) Se puede comprobar cómo el procedimiento para calcular las derivadas es un proceso recurrente, con lo que permite aproximar hasta orden n de forma sencilla. El siguiente paso es la elección de los N nodos, x1, x2, …, xN, en el interior del dominio [a, b], que resulta fundamental para que la función f(x) sea una buena aproximación a la solución del problema en todo el intervalo. Si F(x) es la solución, se puede demostrar por teoría de que existe un punto c en el dominio, tal que se cumple para cualquier punto de dicho dominio que el error entre la aproximación y la solución viene dado por: 42 42 𝑒(𝑥)=𝑓(𝑥)−𝐹(𝑥)=𝐴(𝑥) 𝑁! ·𝑓(𝑁)(𝑐) 𝐴(𝑥)=(𝑥−𝑥1)(𝑥−𝑥2)…(𝑥−𝑥𝑁) (57) Por tanto, conociendo cotas de la derivada enésima de f y para A(x) se podría estimar el error. En el caso de A(x), al tratarse de un polinomio, su valor absoluto tiene un máximo en algún punto del dominio, que dependerá de la elección de los nodos. Resulta ahora evidente que la elección de unos nodos tales que minimicen el valor del máximo de |A(x)| implica la minimización del error cometido, es por ello que se eligen los nodos de Chebyshev, que como se comentó anteriormente, minimizan la norma de A(x). ‖𝐴‖=max 𝑎≤𝑥≤𝑏|𝐴(𝑥)| (58) Método de colocación para problemas bidimensionales con polinomios interpolantes de Lagrange y nodos de Chebyshev Ahora se va un paso más alla, tratando el método de colocación de forma análoga al punto anterior, pero con la diferencia de que se aplica sobre un dominio bidimensional, que será lo que se utilizará en lo que sigue de estudio. Cabe destacar que, en este caso, los nodos elegidos en ambas direcciones son los de Chebyshev. En cambio, la discretización del dominio en los códigos elaborados en este proyecto mezcla nodos de Chebyshev en z y nodos basados en diferencias finitas en x. Por tanto, este apartado continúa siendo una presentación para el lector de las herramientas empleadas, ya que para la resolución del problema se ha tomado un método mixto entre el método de colocación bidimensional con nodos de Chebyshev y el método de colocación bidimensional con diferencias finitas, que se tratará posteriormente. Sea Φ=Φ (𝑥,𝑧) una función definida en un dominio rectangular 0≤𝑥≤𝐿 y 0≤𝑧≤𝐻 en el que se han elegido Nx · Nz nodos { (xm, zn) : m=1 … Nx , n = 1 … Nz }. Para cualquier x dado se puede realizar una interpolación de Lagrange para su dependencia con la coordenada z de la siguiente forma: Φ(𝑥,𝑧)=∑𝐿𝑧,𝑛(𝑧)·Φ(𝑥,𝑧𝑛) 𝑁𝑧 𝑛=1 (60) Si esto se repite para cada una de las funciones de la coordenada x, Φ(𝑥,𝑧𝑛), que aparece en el sumatorio anterior, se tiene la expresión de la fórmula de interpolación de Lagrange para dos variables: Φ(𝑥,𝑧)=∑𝐿𝑧,𝑛(𝑧)·[ 𝑁𝑧 𝑛=1 ∑𝐿𝑥,𝑚(𝑥)·Φ(𝑥𝑚,𝑧𝑛) 𝑁𝑥 𝑚=1 ] (61) Siendo Φ(𝑥𝑚,𝑧𝑛) el valor que toma la función en el nodo (xm,zn), y 𝐿𝑧,𝑛, 𝐿𝑥,𝑚 los polinomios interpolantes de Lagrange para cada variable, que vienen dados por: 𝐿𝑥,𝑚(𝑥)=∏𝑥−𝑥𝑗 𝑥𝑚−𝑥𝑗 𝐿𝑧,𝑛(𝑧)=∏𝑧−𝑧𝑗 𝑧𝑛−𝑧𝑗 𝑁𝑧 𝑗=1 𝑗≠𝑛 𝑁𝑥 𝑗=1 𝑗≠𝑚 (62) Continuando con el cálculo de las derivadas parciales de Φ(𝑥,𝑧), las derivadas parciales con respecto a x y con respecto a z son: Φ𝑥(𝑥,𝑧)≡𝜕Φ(𝑥,𝑧) 𝜕𝑥 =∑∑𝐿𝑥,𝑚 ′(𝑥)·𝐿𝑧,𝑛(𝑧)·Φ(𝑥𝑚,𝑧𝑛) 𝑁𝑧 𝑛=1 𝑁𝑥 𝑚=1 (63) Φ𝑧(𝑥,𝑧)≡𝜕Φ(𝑥,𝑧) 𝜕𝑧 =∑∑𝐿𝑥,𝑚(𝑥)·𝐿𝑧,𝑛 ′(𝑧)·Φ(𝑥𝑚,𝑧𝑛) 𝑁𝑧 𝑛=1 𝑁𝑥 𝑚=1 (64) Llegados a este punto, se eligen los nodos de Chebyshev tanto en dirección x como en z, que son los puntos (xi, zj) del dominio dados por las expresiones: 𝑥𝑖=𝑥min +𝑥𝑚𝑎𝑥 2+𝑥𝑚𝑖𝑛−𝑥𝑚𝑎𝑥 2·cos[𝜋 (𝑖−1) 𝑁𝑥−1] 𝑑𝑜𝑛𝑑𝑒 𝑖=1…𝑁𝑥 (65) 𝑧𝑗=𝑧𝑚𝑖𝑛+𝑧𝑚𝑎𝑥 2+𝑧𝑚𝑖𝑛−𝑧𝑚𝑎𝑥 2·cos[𝜋 (𝑗−1) 𝑁𝑧−1] 𝑑𝑜𝑛𝑑𝑒 𝑗=1…𝑁𝑧 (66) Tras mallar, se enumerará cada nodo con un único índice, lo que que permitir reescribir vectorialmente el problema, para ello se usa una asignación biunívoca. (𝑥𝑖,𝑧𝑗)↔𝐼=(𝑖−1)·𝑁𝑧+𝑗 (67) Se puede observar un ejemplo para Nx = 4 y Nz = 3 en la imagen adjunta a continuación. Ilustración 29 – Ejemplo de numeración de nodos con un índice [18]. Al considerar un único índice se puede llevar a cabo el producto entre los polinomios interpolantes y sus derivadas como elementos de matrices Dx y Dz, definidas por: 𝐷𝑥(𝐼,𝐾)=𝐿𝑥,𝑚 ′(𝑥𝑖)𝐿𝑧,𝑛(𝑧𝑗) (68) 𝐷𝑧(𝐼,𝐾)=𝐿𝑥.𝑚(𝑥𝑖)𝐿𝑧,𝑛 ′(𝑧𝑗) (69) 𝐼,𝐾=1…(𝑁𝑥·𝑁𝑧 = 𝑁𝑡) (70) 𝐼=(𝑖−1)·𝑁𝑧+𝑗 (71) 𝐾=(𝑚−1)·𝑁𝑧+𝑛 (72) Con este nuevo producto de matrices la derivada parcial respecto de x y de z pueden reescribirse de una forma 44 44 más simplificada: 𝜕Φ 𝜕𝑥|𝐼=∑𝐷𝑥(𝐼,𝐾)·Φ𝐾 𝑁𝑡 𝐾=1 (73) 𝜕Φ 𝜕𝑧|𝐼=∑𝐷𝑧(𝐼,𝐾)·Φ𝐾 𝑁𝑡 𝐾=1 (74) Donde Φ𝐾 es el valor de la función en los nodos. Al igual que en el caso unidimensional, se puede hacer uso de vectores para obtener una notación más sencilla y manejable. Φ 󰇍 󰇍 󰇍 ≡[Φ1 Φ2 … Φ𝑁𝑡]; Φ 󰇍 󰇍 󰇍 𝑥≡ [ Φ𝑥1 ′ Φ𝑥2 ′ … Φ𝑥𝑁𝑡 ′ ] ; Φ 󰇍 󰇍 󰇍 𝑧≡ [ Φ𝑧1 ′ Φ𝑧2 ′ … Φ𝑧𝑁𝑡 ′ ] (75) Φ𝑥 󰇍 󰇍 󰇍 󰇍 󰇍 =𝐷𝑥·Φ 󰇍 󰇍 󰇍 Φ𝑧 󰇍 󰇍 󰇍 󰇍 󰇍 =𝐷𝑧·Φ 󰇍 󰇍 󰇍 (76) Para derivadas de órdenes superiores, de forma análoga al problema de una varible, se tiene que: 𝜕2Φ 𝜕𝑥2|𝐼=∑𝐷𝑥(𝐼,𝐾)·Φ𝑥,𝐾 𝑁𝑡 𝐾=1 ≈∑(∑𝐷𝑥(𝐼,𝐾)·Φ𝐾 𝑁𝑡 𝐾=1 )·𝐷𝑥(𝐼,𝐾), 𝐼=1…𝑁𝑡 𝑁𝑡 𝐾=1 (77) Esa misma expresión, se puede conseguir en término de vectores Φ 󰇍 󰇍 󰇍 y Φ𝑥 󰇍 󰇍 󰇍 󰇍 󰇍 y utilizando la matrix 𝐷𝑥, lo que da como resultado la expresión: Φ𝑥𝑥 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 =𝐷𝑥·Φ𝑥 󰇍 󰇍 󰇍 󰇍 󰇍 =𝐷𝑥2Φ 󰇍 󰇍 󰇍 (78) Siendo las componentes del vector Φ𝑥𝑥 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 los valores aproximados de Φ𝑥𝑥 en los puntos del mallado sobre el dominio. Un proceso idéntico se sigue para las derivadas de orden superior en z, y para la derivada mixta, el resultado es: Φ𝑧𝑧 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 =𝐷𝑧·Φ𝑧 󰇍 󰇍 󰇍 󰇍 󰇍 =𝐷𝑧2Φ 󰇍 󰇍 󰇍 (79) Φ𝑥𝑧 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 =𝐷𝑥·𝐷𝑧·Φ 󰇍 󰇍 󰇍 =𝐷𝑧·𝐷𝑥·Φ 󰇍 󰇍 󰇍 =Φ𝑧𝑥 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 (80) Método de colocación para problemas bidimensionales con polinomios interpolantes de Lagrange y diferencias finitas El procedimiento es análogo al tratado en los casos anteriores. Al igual que con el caso con nodos de Chebyshev, las matrices derivadas se calculan utilizando polinomios interpolantes de Lagrage según la implementación en los códigos adjuntos al final de la sección 2.7. Cabe mencionar que se entiende como “diferencias finitas” a la consideración de un entorno para cada punto de la malla que, a diferencia del caso con Chebyshev, contenga un número de puntos menor al de la malla total, en este caso un entorno de siete puntos, que permite el cálculo hasta la derivada cuarta con una precisión de ℎ3. El resultado son matrices derivadas dispersas. Ello hace que esta metodología sea óptima para el caso de mallas con un gran número de puntos. Así, para el código complejo, donde la precisión mejora, se utiliza este método en ambas direcciones. 2.4 Particularización al problema eléctrico Tal y como se trató durante el apartado 2.1 – Formulación del problema completo, la resolución del problema eléctrico está desacoplada de la resolución del problema fluido, lo cual permite obtener la fuerza electrohidrodinámica e incluirla posteriormente en las ecuaciones de Navier-Stokes. La obtención de la fuerza eléctrica pasa por calcular el potencial asociado a campos eléctricos externos y la densidad de carga sobre el electrodo embebido, según las expresiones (15) y (25). Para su resolución se aplicará el método de colocación. Pueden consultarse más detalles sobre los métodos utilizados para la resolución numérica del problema en [43]. Cálculo del potencial eléctrico externo Se presenta un problema gobernado por una ecuación de Poisson, donde es necesario calcular el potencial a lo largo de dos regiones, la fluida y la dieléctrica, que presentan distinta permitividad relativa. Una ecuación de Poisson viene dada por la fórmula: 𝜕2Φ 𝜕𝑥2+𝜕2Φ 𝜕𝑧2=𝑓(𝑥,𝑧) (81) Para un dominio rectangular de altura H y longitud L y siendo el término fuente, f (x, z), conocido. La ecuación a resolver para el potencial es la siguiente, ∇(𝜀𝑟·∇𝜙)=0 (82) Expresión que se puede reescribir como: 𝜕2𝜙 𝜕𝑥2+𝜕2𝜙 𝜕𝑧2=0 (83) El dominio es completo, un dominio rectangular que acopla la región fluida, comprendida entre 𝑥𝑚𝑖𝑛,𝑥𝑚𝑎𝑥,𝑧𝑎,𝑚𝑖𝑛 y 𝑧𝑎,𝑚𝑎𝑥, y la región dieléctrica, delimitada por 𝑥𝑚𝑖𝑛,𝑥𝑚𝑎𝑥,𝑧𝑏,𝑚𝑖𝑛 y 𝑧𝑏,𝑚𝑎𝑥. Las condiciones de contorno necesarias para la resolución del sistema aparecen reflejadas en la siguiente imagen: Ilustración 30 – Condiciones de contorno para el cálculo del potencial. Las condiciones de contorno en la entrefase, relacionaran dos puntos infinitamente próximos, cada uno a un lado de la superficie de separación entre los dos medios, medios que a su vez son eléctricamente distintos ya que no están caracterizados por la misma conductividad, ni por la misma permitividad o permeabilidad magnética. Así, para cualquier problema de campo electromagnético, deben satisfacerse una serie de Aire: ∇ (𝜀𝑟1·∇𝜙 )=0 𝜀𝑟1=1 Contorno exterior: ∂ϕ ∂n=0 𝜙=𝜙(𝑡) 𝜙=0 𝜀𝑟=𝑚𝑒𝑑𝑖𝑎(𝜀𝑟1,𝜀𝑟2) Dieléctrico (Kapton): ∇ (𝜀𝑟2·∇𝜙 )=0 𝜀𝑟2=2.7 46 46 condiciones de contorno sobre la entrefase: 𝑛𝑎𝑏 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 ·(𝜀𝑏·𝐸 󰇍 𝑏−𝜀𝑎·𝐸𝑎 󰇍 󰇍 󰇍 󰇍 )=𝜌𝑠 (84) 𝑛𝑎𝑏 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 ·(𝐵 󰇍 𝑏−𝐵𝑎 󰇍 󰇍 󰇍 󰇍 )=0 (85) 𝑛𝑎𝑏 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 ·(𝐽𝑏 󰇍 󰇍 󰇍 −𝐽𝑎 󰇍 󰇍 󰇍 )=−𝜕𝜌𝑐 𝜕𝑡 (86) 𝑛𝑎𝑏 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 ×(𝐻𝑏 󰇍 󰇍 󰇍 󰇍 󰇍 −𝐻𝑎 󰇍 󰇍 󰇍 󰇍 󰇍 )=𝐽 𝑆 (87) 𝑛𝑎𝑏 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 ×(𝐸𝑏 󰇍 󰇍 󰇍 󰇍 −𝐸𝑎 󰇍 󰇍 󰇍 󰇍 )=0 (88) Donde 𝑛𝑎𝑏 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 es el vector normal a la entrefase que apunta desde el dominio a al b, 𝜀 es la permitividad dieléctrica del medio, 𝐸 󰇍 el campo eléctrico, 𝜌𝑠 es la densidad de carga superficial en la entrefase, 𝐵 󰇍 es el campo magnético, 𝐽 es la densidad de corriente, 𝐽𝑠 󰇍 󰇍 la densidad de corriente en la entrefase y 𝐻 󰇍 󰇍 es la intensidad de campo magnético. Si tenemos en cuenta que la variación temporal de las fuerzas magnéticas es despreciable en el caso de los actuadores de plasma, las condiciones de contorno del problema se relajan, quedando resumidas en: 𝑛𝑎𝑏 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 ·(𝜀𝑏·𝐸 󰇍 𝑏−𝜀𝑎·𝐸𝑎 󰇍 󰇍 󰇍 󰇍 )=𝜌𝑠 (89) 𝑛𝑎𝑏 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 ×(𝐸𝑏 󰇍 󰇍 󰇍 󰇍 −𝐸𝑎 󰇍 󰇍 󰇍 󰇍 )=0 (90) La primera de ellas impone que, si en la frontera entre los dos medios existe densidad de carga superficial distinta de cero, entonces la componente normal del desplazamiento eléctrico (𝜀·𝐸) sufre un salto no continuo igual a la densidad de carga superficial en la entrefase. Si no hay densidad de carga superficial en la entrefase, entonces el desplazamiento eléctrico sí es contínuo. En este caso, la densidad de carga superficial es nula, lo cual se podría comprobar porque sobre un elemento de volumen en la entrefase (𝑑𝑉=𝑑𝐴·𝑑𝑧=𝑑𝑥·𝑑𝑦·𝑑𝑧) la carga total vendría dada por 𝜌𝑐·𝑑𝑉, siendo la densidad de carga finita. La carga por unidad de superficie vendría dada por 𝜌𝑐·𝑑𝑉 𝑑𝐴=𝜌𝑐·𝑑𝑧, donde 𝑑𝑧→0, lo que implicaría que la densidad de carga superficial es nula. La segunda condición garantiza que la componente tangencial del campo eléctrico no varía a un lado y a otro de la entrefase. Por consiguiente, 𝜕𝜙 𝜕𝑛(𝑥𝑚𝑖𝑛,𝑧)=0, 𝜕𝜙 𝜕𝑛(𝑥𝑚𝑎𝑥,𝑧)=0, 𝜕𝜙 𝜕𝑛(𝑥,𝑧𝑎,𝑚𝑎𝑥)=0, 𝜕𝜙 𝜕𝑛(𝑥,𝑧𝑏,𝑚𝑖𝑛)=0, 𝜙(𝑥𝑒𝑙𝑒𝑐𝑡 𝑜𝑐𝑢𝑙𝑡𝑜,𝑧𝑏,𝑚𝑖𝑛)=0 𝜙(𝑥𝑒𝑙𝑒𝑐𝑡 𝑣𝑖𝑠𝑡𝑜,𝑧𝑎,𝑚𝑖𝑛)=𝜙(𝑥𝑒𝑙𝑒𝑐𝑡 𝑣𝑖𝑠𝑡𝑜,𝑧𝑏,𝑚𝑎𝑥)=𝜙𝑚𝑎𝑥 𝐸𝑥𝑎(𝑥,𝑧𝑎,𝑚𝑖𝑛)=𝐸𝑥𝑏(𝑥,𝑧𝑏,𝑚𝑎𝑥)≡−𝜕𝜙 𝜕𝑥(𝑥,𝑧𝑎,𝑚𝑖𝑛)=−𝜕𝜙 𝜕𝑥(𝑥,𝑧𝑏,𝑚𝑎𝑥) 𝜀𝑎𝐸𝑧𝑎(𝑥,𝑧𝑎,𝑚𝑖𝑛)=𝜀𝑏𝐸𝑧𝑏(𝑥,𝑧𝑏,𝑚𝑎𝑥)≡−𝜀𝑎𝜕𝜙 𝜕𝑧(𝑥,𝑧𝑎,𝑚𝑖𝑛)=−𝜀𝑏𝜕𝜙 𝜕𝑧(𝑥,𝑧𝑏,𝑚𝑎𝑥) (91) Resulta más conveniente resolver el problema de manera adimensional, para ello se definen las siguientes variables adimensionales: 𝑥=𝑥𝐿, 𝑧=𝑧𝐿, 𝜏=𝑡𝑡0, 𝜙∗=𝜙 𝜙𝑚𝑎𝑥·𝑓(𝑡), 𝐸∗ 󰇍 󰇍 󰇍 󰇍 =𝐿 𝜙𝑚𝑎𝑥𝑓(𝑡)·𝐸 󰇍 (92) Donde 𝐿 es la longitud característica del problema, basada en [19] y que se corresponde con el gap entre los dos electrodos, en este caso 𝐿=0.001 𝑚, y 𝜙𝑚𝑎𝑥 es la amplitud del voltaje aplicado sobre los electrodos, que se toma igual a 5 kV, como en la publicación de Suzen et al. [15]. El problema adimensional resultante es el siguiente: 𝜕2𝜙∗ 𝜕𝑥2+𝜕2𝜙∗ 𝜕𝑧2=0 (93) 𝜕𝜙∗ 𝜕𝑛(𝑥𝑚𝑖𝑛,𝑧)=0, 𝜕𝜙∗ 𝜕𝑛(𝑥𝑚𝑎𝑥,𝑧)=0, 𝜕𝜙∗ 𝜕𝑛(𝑥,𝑧𝑎,𝑚𝑎𝑥)=0, 𝜕𝜙∗ 𝜕𝑛(𝑥,𝑧𝑏,𝑚𝑖𝑛)=0, 𝜙∗(𝑥𝑒𝑙𝑒𝑐𝑡 𝑜𝑐𝑢𝑙𝑡𝑜,𝑧𝑏,𝑚𝑖𝑛)=0 𝜙∗(𝑥𝑒𝑙𝑒𝑐𝑡 𝑣𝑖𝑠𝑡𝑜,𝑧𝑎,𝑚𝑖𝑛)=𝜙∗(𝑥𝑒𝑙𝑒𝑐𝑡 𝑣𝑖𝑠𝑡𝑜,𝑧𝑏,𝑚𝑎𝑥)=1 𝐸𝑥𝑎 ∗(𝑥,𝑧𝑎,𝑚𝑖𝑛)=𝐸𝑥𝑏 ∗(𝑥,𝑧𝑏,𝑚𝑎𝑥)≡−𝜕𝜙∗ 𝜕𝑥(𝑥,𝑧𝑎,𝑚𝑖𝑛)=−𝜕𝜙∗ 𝜕𝑥(𝑥,𝑧𝑏,𝑚𝑎𝑥) 𝜀𝑎𝐸𝑧𝑎 ∗(𝑥,𝑧𝑎,𝑚𝑖𝑛)=𝜀𝑏𝐸𝑧𝑏 ∗(𝑥,𝑧𝑏,𝑚𝑎𝑥)≡−𝜀𝑎𝜕𝜙∗ 𝜕𝑧(𝑥,𝑧𝑎,𝑚𝑖𝑛)=−𝜀𝑏𝜕𝜙∗ 𝜕𝑧(𝑥,𝑧𝑏,𝑚𝑎𝑥) (94) El dominio se divide en dos regiones, cada una de ellas malladas con 𝑁𝑡𝑎=𝑁𝑥·𝑁𝑧𝑎 y 𝑁𝑡𝑏=𝑁𝑥·𝑁𝑧𝑏 puntos respectivamente, por tanto, el número de nodos totales es 𝑁𝑡=𝑁𝑡𝑎+𝑁𝑡𝑏. Tanto la región fluida como la dieléctrica presentan mallados con nodos de Chebyshev en la dirección z en el caso del código sencillo, y nodos basados en diferencias finitas para el código más complejo. Además, en ambos casos se usan nodos basados en diferencias finitas en x. De nuevo, enumerando cada uno de los nodos con un solo índice en lugar de dos se puede expresar el producto de los polinomios interpoladores y de sus derivadas como elementos de matrices 𝐷𝑥,𝐷𝑧 y ello permite reescribir el problema de forma vectorial. (𝐷𝑥2+𝐷𝑧2) 𝜙∗ 󰇍 󰇍 󰇍 󰇍 ≡𝐷𝐿 𝜙∗ 󰇍 󰇍 󰇍 󰇍 =0 󰇍 (95) Donde el vector 𝜙 󰇍 contiene el valor del potencial en cada uno de los nodos de la malla. Llegados a este punto, resulta necesario acoplar las dos regiones para resolver el problema completo con un sólo sistema, que a su vez tenga en cuenta los aspectos asociados a la región fluida, los asociados a la región dieléctrica y la interacción entre ambos medios a través de la entrefase común. Mediante un tratamiento sistemático, que se aplicará también en el cálculo de la densidad de carga y que puede generalizarse para acoplar tantas regiones como se desee, se confecciona el sistema de ecuaciones que permite resolver el problema completo. La metodología seguida en este estudio no está presente en ninguna de las publicaciones consultadas, y es por ello que quede aquí recogida. Sea un sistema de la forma: 𝐴·𝑥 =𝑏 󰇍 (96) En primer lugar, se construye una matriz del sistema, 𝐴, que a su vez está formada por varias submatrices: 𝐴=[𝐴11 𝐴12 ⋯ 𝐴1𝑁 𝐴21 𝐴22 𝐴23 ⋮ ⋮ 𝐴32 ⋱ 𝐴𝑁−1,𝑁 𝐴𝑁1 ⋯ 𝐴𝑁,𝑁−1 𝐴𝑁𝑁 ] (97) Las matrices de la diagonal están asociadas directamente con cada uno de los medios por separado, el resto de matrices recogen el acoplamiento y la interacción entre las distintas regiones, dichas relaciones se incluyen a posteriori, al imponer las condiciones de contorno sobre el sistema. El siguiente paso es construir el término independiente del sistema, para ello se crea un vector columna que agrupe en su interior los términos independientes de los subsistemas correspondientes a cada dominio, es 54 54 Implementación del problema de espesor en el cálculo de la densidad de carga El código más complejo incorpora el espesor del electrodo visto, así se incluye la influencia de este sobre los resultados obtenidos al calcular la densidad de carga. Ilustración 33 – Nodos que contienen al electrodo visto. Para tener en cuenta el espesor de la lámina conductora se localizan los nodos que la contienen en función de la altura del electrodo. Una vez determinados qué nodos están asociados al electrodo visto se imponen las siguientes condiciones de contorno sobre ellos. Ilustración 34 – Condiciones de contorno adimensionales impuestas sobre el electrodo visto. Cálculo de la fuerza electrohidrodinámica Una vez que el potencial eléctrico y la densidad de carga son conocidas las componentes en x y z de la fuerza generada por el actuador se calculan según las siguientes expresiones: 𝐹𝑏𝑥 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 =𝜌𝑐 󰇍 󰇍 󰇍 󰇍 ·𝐷𝑥·𝜙 󰇍 =𝜌𝑐 󰇍 󰇍 󰇍 󰇍 ·𝐸𝑥 󰇍 󰇍 󰇍 󰇍 (125) 𝐹𝑏𝑧 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 =𝜌𝑐 󰇍 󰇍 󰇍 󰇍 ·𝐷𝑧·𝜙 󰇍 =𝜌𝑐 󰇍 󰇍 󰇍 󰇍 ·𝐸𝑧 󰇍 󰇍 󰇍 󰇍 (126) Aire Dieléctrico (Kapton) 𝜌𝑐∗=0 𝜌𝑐∗=0 𝜌𝑐∗=0 𝜌𝑐∗=0 2.5 Particularización al problema fluido Tras incorporar la fuerza electrohidrodinámica al problema fluido, las ecuaciones de Navier-Stokes para el caso incompresible bidimensional se resolverán mediante la función de corriente, que como se trató en el apartado 2.1 – Formulación del problema completo, permiten pasar a una única ecuación con una incógnita. En este proyecto queda recogido un breve resumen sobre el proceso que se llevó a cabo para la resolución del problema hidrodinámico, sin embargo, se puede encontrar un desarrollo mucho más extenso y detallado en el trabajo de fin de grado de Dugo, P. [9]. Cálculo del campo de velocidades La ecuación del problema fluido es de tipo biharmónica como se puede comprobar a continuación: 𝜕 𝜕𝜏(∇2Ψ)+𝜕Ψ 𝜕𝑧1 𝜕𝑥 [∇2Ψ]−𝜕Ψ 𝜕𝑥1 𝜕𝑧[∇2Ψ]−1 𝑅𝑒∇4Ψ=𝐷𝑐(𝜕(𝜌𝑐∗𝐸𝑥∗) 𝜕𝑧 −𝜕(𝜌𝑐∗𝐸𝑧∗) 𝜕𝑥 )𝑓(𝜏)2 (127) El dominio sobre el que actúa es una región rectangular de altura comprendida entre 𝑧𝑎,𝑚𝑖𝑛 y 𝑧𝑎,𝑚𝑎𝑥, y longitud comprendida entre 𝑥𝑎,𝑚𝑖𝑛 y 𝑥𝑎,𝑚𝑎𝑥. Para la resolución numérica de la ecuación es necesario discretizar la región fluida mediante un mallado con 𝑁𝑡𝑎=𝑁𝑥𝑎·𝑁𝑧𝑎 nodos. Además, la expresión se resuelve sujeta a una serie de condiciones de contorno recogidas a continuación: Ilustración 35 – Condiciones de contorno para resolver el problema fluido. Ψ(𝑥𝑚𝑖𝑛,𝑧)=0, Ψ(𝑥𝑚𝑎𝑥,𝑧)=0, Ψ(𝑥,𝑧𝑎,𝑚𝑖𝑛)=0 , Ψ(𝑥,𝑧𝑎,𝑚𝑎𝑥)=0 (128) ∂Ψ ∂x(𝑥𝑚𝑖𝑛,𝑧)=0, ∂Ψ ∂x(𝑥𝑚𝑎𝑥,𝑧)=0, ∂Ψ ∂z(𝑥,𝑧𝑎,𝑚𝑖𝑛)=0 , ∂Ψ ∂z(𝑥,𝑧𝑎,𝑚𝑎𝑥)=0 (129) En la ecuación está presente la variable tiempo, por lo que resulta necesario discretizar las derivadas temporales para poder continuar aplicando el método de colocación: 𝜕 𝜕𝜏(∇2Ψ)=∇2Ψ(𝑥,𝑧,𝜏𝑛)−∇2Ψ(𝑥,𝑧,𝜏𝑛−1) τn−𝜏𝑛−1 =∇2Ψ𝑛−∇2Ψ𝑛−1 Δ𝜏 (130) Donde Ψ𝑛 y Ψ𝑛−1 son la función de corriente en un instante y en el instante anterior respectivamente. La ecuación se podría reescribir como: (∇2Ψn−∇2Ψn−1)+(𝜕Ψ 𝜕𝑧1 𝜕𝑥 [∇2Ψ]−𝜕Ψ 𝜕𝑥1 𝜕𝑧[∇2Ψ]−1 𝑅𝑒∇4Ψ)·Δτ=𝐷𝑐(𝜕(𝜌𝑐∗𝐸𝑥∗) 𝜕𝑧 −𝜕(𝜌𝑐∗𝐸𝑧∗) 𝜕𝑥 )𝑓(𝜏)2·Δ𝜏 (131) Aire Ψ=𝜕Ψ 𝜕𝑛=0 Dieléctrico (Kapton) Ψ=𝜕Ψ 𝜕𝑛=0 56 56 Enumerando cada uno de los nodos del dominio mediante un sólo índice con la expresión (58) se puede expresar el producto entre los polinomios interpolantes del método de colocación y sus derivadas como elementos de las matrices 𝐷𝑥𝑎,𝐷𝑧𝑎. Esto último permite aproximar el valor de la derivada parcial respecto x o z de una magnitud como el producto de esta por una de las matrices anteriores, y expresar la ecuación del problema de forma vectorial. La expresión del operador biharmónico de la función de corriente se puede calcular como: ∇4 Ψ 󰇍 󰇍 󰇍 󰇍 =(𝐷𝑥𝑎 2+𝐷𝑧𝑎 2)·(𝐷𝑥𝑎 2+𝐷𝑧𝑎 2)·Ψ 󰇍 󰇍 󰇍 =(𝐷𝑥𝑎 4+𝐷𝑧𝑎 4+2·𝐷𝑥𝑎 2𝐷𝑧𝑎 2)·Ψ 󰇍 󰇍 󰇍 =𝐷𝐿2·Ψ 󰇍 󰇍 󰇍 (132) Por consiguiente, la ecuación completa expresada de forma vectorial sería: 𝐷𝐿(Ψn 󰇍 󰇍 󰇍 󰇍 󰇍 −Ψ 󰇍 󰇍 󰇍 𝑛−1)+(𝐷𝑧𝑎Ψ 󰇍 󰇍 󰇍 𝑛·𝐷𝑥𝑎𝐷𝐿Ψ 󰇍 󰇍 󰇍 𝑛−𝐷𝑥𝑎Ψ 󰇍 󰇍 󰇍 𝑛·𝐷𝑧𝑎𝐷𝐿Ψ 󰇍 󰇍 󰇍 𝑛)·Δ𝜏−1 𝑅𝑒𝐷𝐿2Ψ 󰇍 󰇍 󰇍 𝑛·Δ𝜏 =𝐷𝑐·(𝐷𝑧𝑎·(𝜌𝑐∗𝐸𝑥∗)−𝐷𝑥𝑎·(𝜌𝑐∗𝐸𝑧∗))·𝑓(𝜏)2·Δ𝜏=𝐹 (133) Donde los vectores Ψ 󰇍 󰇍 󰇍 y 𝐹 contienen los valores de Ψ(𝑥,𝑧,𝜏) y 𝐹(𝑥,𝑧,𝜏) en los nodos de la malla. Siendo 𝐷𝐿2 la matriz del operador biharmónico y 𝐹, el término fuente, que es conocido. Además, parte de la ecuación aparece recuadrada, es el término no lineal (TNL), y para calcularlo se aproxima mediante el mismo término en el instante anterior, es decir: 𝑇𝑁𝐿𝑛≈𝑇𝑁𝐿𝑛−1 (134) Expresando el sistema como: 𝐴·Ψ 󰇍 󰇍 󰇍 𝑛=𝑏 󰇍 (135) La matriz del sistema sería: 𝐴=𝐷𝐿−1 𝑅𝑒·𝐷𝐿2·Δ𝜏 (136) Y el término independiente: 𝑏 󰇍 =𝐹 − 𝑇𝑁𝐿 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 𝑛−1·Δ𝜏+𝐷𝐿·Ψ 󰇍 󰇍 󰇍 𝑛−1 (137) Por último, antes de resolver el sistema se deben incluir las condiciones de contorno. Para las condiciones correspondientes a los valores de la función Ψ sobre las paredes del dominio se procede de la misma forma que con el problema eléctrico. Se modifica la matriz del sistema, en este caso 𝐴, y el vector de términos independientes, 𝑏 󰇍 . Si el valor del índice I coincide con el de un nodo situado sobre alguna pared los elementos de la fila I de la matriz 𝐴 se anulan, excepto en la diagonal, que toma valor unidad. Además, el elemento 𝑏 󰇍 (𝐼) toma el valor correspondiente a la condición de contorno impuesta, en este caso 𝑏 󰇍 (𝐼)=0, tanto si el nodo pertenece a una pared vertical como a una pared horizontal. Para el caso de las condiciones de contorno asociadas a derivadas parciales de la función Ψ, estas se imponen sobre los nodos del subcontorno adyacente a las paredes del dominio.  Si I coincide con un nodo perteneciente al subcontorno adyacente a la pared vertical izquierda, se sustituye la fila I de 𝐴 por la fila de la matriz 𝐷𝑥𝑎 correspondiente al nodo K=I-Nz de la pared vertical izquierda, es decir, se hace 𝐴(𝐼,1:𝑁𝑡)=𝐷𝑥𝑎(𝐼−𝑁𝑧,1:𝑁𝑡). De forma análoga, para el margen vertical derecho, el cambio sería 𝐴(𝐼,1:𝑁𝑡)=𝐷𝑥𝑎(𝐼+𝑁𝑧,1:𝑁𝑡). Además, se impone el valor del término independiente en la componente I según la condición de contorno elegida, en este caso, 𝑏 󰇍 (𝐼)=0.  Si I coincide con un nodo perteneciente al subcontorno adyacente a la pared horizontal inferior, se sustituye la fila I de 𝐴 por la fila de la matriz 𝐷𝑧𝑎 correspondiente al nodo K=I-1 de la pared horizontal inferior, es decir, se hace 𝐴(𝐼,1:𝑁𝑡)=𝐷𝑧𝑎(𝐼−1,1:𝑁𝑡). Se procede de forma equivalente para la pared horizontal superior, con el cambio 𝐴(𝐼,1:𝑁𝑡)=𝐷𝑧𝑎(𝐼+1,1:𝑁𝑡). Finalmente, se impone el valor de la componente I en el término independiente en función de la condición de contorno impuesta, en este caso, 𝑏 󰇍 (𝐼)=0. Ilustración 36 – Contorno y subcontorno del dominio. 2.6 Problema con condiciones de contorno periódicas El modelado del actuador de plasma se ha llevado a cabo con condiciones de contorno estacionarias, sin embargo, resulta interesante particularizarlo al caso de condiciones periódicas sobre paredes verticales del dominio. La imposición de condiciones periódicas en las paredes verticales hace que el sistema que se está estudiando esté flanqueado a cada lado por infinitos sistemas idénticos. Esta característica hace más flexible al modelo de actuador elaborado, permitiendo estudiar el efecto de acoplar infinitos dispositivos idénticos en línea y modificar la distancia entre ellos con sólo variar la longitud del dominio rectangular considerado. Así, se puede estudiar cómo es la interacción entre actuadores suficientemente cercanos. Es un problema interesante desde el punto de vista eléctrico, permitiendo estudiar cómo varía el potencial y la densidad de carga ante la presencia de otros actuadores, sin embargo, sobre el problema fluido estas condiciones no presentan gran interés puesto que al igualar las condiciones de entrada en la pared izquierda y en la derecha no se puede estudiar cómo es la evolución del fluido al pasar por actuadores sucesivos. Como se comentó anteriormente, este tipo de condiciones se caracteriza por igualar en las dos paredes verticales del dominio la magnitud que se está estudiando y su derivada respecto a la dirección normal. A continuación, se desarrolla la forma en la que se han implantado sobre el problema eléctrico:  Cálculo de la densidad de carga: 𝜌𝑐∗(𝑥𝑚𝑖𝑛,𝑧)=0, 𝜌𝑐∗(𝑥𝑚𝑎𝑥,𝑧)=0→ 𝜌𝑐∗(𝑥𝑚𝑖𝑛,𝑧)=𝜌𝑐∗(𝑥𝑚𝑎𝑥,𝑧),𝜕𝜌𝑐∗ 𝜕𝑛(𝑥𝑚𝑖𝑛,𝑧)=𝜕𝜌𝑐∗ 𝜕𝑛(𝑥𝑚𝑎𝑥,𝑧) (138)  Cálculo del potencial eléctrico debido a cargas externas: 𝜕𝜙∗ 𝜕𝑛(𝑥𝑚𝑖𝑛,𝑧)=0, 𝜕𝜙∗ 𝜕𝑛(𝑥𝑚𝑎𝑥 𝐿,𝑧)=0→ 𝜙∗(𝑥𝑚𝑖𝑛,𝑧)=𝜙∗(𝑥𝑚𝑎𝑥,𝑧),𝜕𝜙∗ 𝜕𝑛(𝑥𝑚𝑖𝑛,𝑧)=𝜕𝜙∗ 𝜕𝑛(𝑥𝑚𝑎𝑥,𝑧) (139) 58 58 La implementación de las condiciones de contorno en el código numérico es la misma para ambos problemas, se imponen tanto sobre la matriz del dominio fluido como sobre la del dominio dieléctrico, aunque aquí se recoge el caso para un subdominio cualquiera. Para ello se define 𝐼𝑖𝑧𝑞 como el índice que coincide con los nodos de la pared vertical izquierda del dominio e 𝐼𝑑𝑒𝑟 como el índice de los nodos coincidentes con la pared vertical derecha: 𝐴(𝐼𝑖𝑧𝑞,:)=0, 𝐴(𝐼𝑖𝑧𝑞,𝐼𝑖𝑧𝑞)=1, 𝐴(𝐼𝑖𝑧𝑞,𝐼𝑑𝑒𝑟)=−1, 𝑏(𝐼𝑖𝑧𝑞,1)=0 𝐴(𝐼𝑑𝑒𝑟,:)=𝐷𝑥(𝐼𝑑𝑒𝑟,:)−𝐷𝑥(𝐼𝑖𝑧𝑞,:), 𝑏(𝐼𝑑𝑒𝑟,1)=0, (140) 2.7 Descripción de los códigos elaborados Este trabajo está basado principalmente en la implementación del problema eléctrico del modelo de actuador de plasma, por ello a continuación se adjuntan los códigos creados para el cálculo de la densidad de carga sobre el dieléctrico y del potencial asociado al campo eléctrico aplicado sobre los electrodos. Aunque aquí se incluye también el código correspondiente al problema fluido, la justificación del mismo queda recogida en el proyecto de Dugo, P. [9]. Definición de variables Es necesario definir una serie de parámetros que se recogen en tres scripts:  Parámetros_mallado: function [Nx,Nza,Nzb,Nta,Ntb,xchad,za_ad,zb_ad,L]=parametros_mallado %PARÁMETROS MALLADO L=0.001; %Longitud característica del problema [m] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% DOMINIO GENERAL %%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % (La X se mantiene igual en los dos dominios) xminad=-0.03/L; %Posición mínima adimensional x del dominio fluido xmaxad=0.06/L; %Posición máxima adimensional x del dominio fluido Nx=350; %Número de puntos en x x1=0.005/L; %Primera división del refinado en x de la malla x2=0.03/L; %Segunda división del refinado en x de la malla N1=50; %Número de puntos en el primer segmento del mallado N2=250; %Número de puntos en el segundo segmento del mallado N3=50; %Número de puntos en el tercero segmento del mallado %Distribución Diferencias finitas de nodos en x xchad(1:N1)=linspace(xminad,x1,N1); xchad((N1):(N2+N1))=linspace(x1,x2,N2+1); xchad((N1+N2):Nx)=linspace(x2,xmaxad,N3+1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% DOMINIO FLUIDO %%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% zaminad=0/L; %Posición mínima adimensional z del dominio fluido zamaxad=0.01/L; %Posición máxima adimensional z del dominio fluido Nza=120; %Número de puntos en z z1=0.002/L; N1=80; N2=40; %Distribución Diferencias finitas de nodos en z za_ad(1:N1)=linspace(zaminad,z1,N1); za_ad(N1:Nza)=linspace(z1,zamaxad,N2+1); %Distribución según x^2 de nodos en z %psi=linspace(0,1,Nza); %za_ad=zamaxad*psi.^2; %Distribución Chebyshev de nodos en z %za_ad(1:Nza)=(zamaxad+zaminad)/2-(zamaxad-zaminad)/2*cos(((1:Nza)-1)*pi/(Nza-1)); Nta=Nx*Nza; %Número total de nodos en el dominio fluido %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% DOMINIO DIELÉCTRICO %%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% zbminad=-0.000127/L; %Posición mínima adimensional z del dominio dieléctrico zbmaxad=0/L; %Posición máxima adimensional z del dominio dieléctrico Nzb=30; %Número de puntos en z %Distribución Diferencias finitas de nodos en z zb_ad(1:Nzb)=linspace(zbminad,zbmaxad,Nzb)'; %Distribución Chebyshev de nodos en z %zb_ad(1:Nzb)=(zbmaxad+zbminad)/2-(zbmaxad-zbminad)/2*cos(((1:Nzb)-1)*pi/(Nzb-1)); Ntb=Nx*Nzb; %Número total de nodos en el dominio fluido Se establecen los límites espaciales del dominio, así como el número de nodos deseados o cómo es la propia distribución de estos en cada una de las regiones. Puede observarse que la distribución en z según Chebyshev aparece comentada, sólo sería necesario descomentarla y comentar la distribución con diferencias finitas para aplicarla. 60 60  Parámetros_eléctricos: function [eps_kapton,eps_aire,ldebye,L,mu,sigma]=parametros_electrico %PARÁMETROS PROBLEMA ELÉCTRICO eps_kapton=2.7; %Permitividad dielectrico (Kapton) eps_aire=1; %Permitividad aire ldebye=0.0001; %Longitud de Debye [m]; La distribución de densidad de carga es bastante sensible con esta longitud L=0.001; %Longitud característica [m] rhoc_max=0.0008; %Densidad de carga máxima, [C/m^3] phi_max=5e3; %Potencial máximo mu=0.0205/L; %Factor posicionamiento distribución Gauss sigma=0.3*(0.02)/L; %Factor escala distribución Gauss, evalua la evolución de la densidad de carga Se introducen los valores de las magnitudes que influyen en el problema eléctrico.  Parámetros_fluidos: function [dtau,u0,Re,Dc]=parametros_fluido L=0.001; %Longitud característica dtau=0.0001; %Paso del tiempo rhoc_max=0.0008; %Densidad de carga máxima, [C/m^3] phi_max=5e3; %Potencial máximo rho=1.225; %Densidad [kg/m^3] mu=1.85e-5; %Viscosidad dinámica u0=sqrt(rhoc_max*phi_max/rho); %Velocidad característica Re=rho*u0*L/mu; %Número Reynolds Dc=rhoc_max*phi_max/u0^2/rho; %Número adimensional fuerza Coulombiana De nuevo se introducen las magnitudes necesarias para el problema, esta vez el fluido. Cálculo de la densidad de carga function DensidadCarga_2Domains close all tic % EXTREMOS ELECTRODOS ini_ev=0.01; fin_ev=0.02; ini_eo=0.0205; fin_eo=0.0305; % CARGA DE PARÁMETROS Y DATOS [Nx,Nza,Nzb,Nta,Ntb,x_ad,za_ad,zb_ad]=parametros_mallado; [eps_kapton,eps_aire,ldebye,L,mu,sigma]=parametros_electrico; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% DOMINIO FLUIDO %%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % MATRICES DERIVADAS [dx, dx2,~,~]=Derivs(x_ad,Nx); [dz,dz2,~,~]=Derivs(za_ad,Nza); Dxa=sparse(kron(dx,speye(Nza))); Dxa2=sparse(kron(dx2,speye(Nza))); Dza=sparse(kron(speye(Nx),dz)); Dza2=sparse(kron(speye(Nx),dz2)); DLa=sparse(Dxa2+Dza2) ; % PERMITIVIDAD RELATIVA DEL MEDIO eps_a=ones(Nta,1); eps_mean=(2*eps_aire*eps_kapton)/(eps_aire+eps_kapton); %Media armónica for i=1:Nx eps_a((i-1)*Nza+1,1)=eps_mean; %permitividad en la entrefase end % MATRIZ CON LA ECUACIÓN A RESOLVER Aa=sparse(eps_a.*DLa-speye(Nta)*(L^2/ldebye^2)); % TÉRMINO INDEPENDIENTE DEL SISTEMA Aa·rhoca=ba ba(1:Nta,1)=0; % CONDICIONES DE CONTORNO % 1) Contornos verticales: rho_c=0; for j=1:Nza I=j; Aa(I,:)=0; Aa(I,I)=1; ba(I,1)=0; I=(Nx-1)*Nza+j; Aa(I,:)=0; Aa(I,I)=1; ba(I,1)=0; end % 2) Contornos horizontal superior: for i=2:(Nx-1) I=(i-1)*Nza+Nza; Aa(I,:)=0; Aa(I,I)=1; ba(I,1)=0; end % PROBLEMA ESPESOR ELECTRODO VISTO %Posición en X de los nodos j=1; for i=1:Nx-1 62 62 if (x_ad(i)-ini_ev/L)*(fin_ev/L-x_ad(i)) >= 0 %Busca el punto de la malla más similar a la posición donde quieres establecer el electrodo (por debajo siempre). ix_elect(j)=i; j=j+1; end end %Determinación de los nodos en Z para dar espesor espesor_e=0.0001/L; %Espesor del electrodo visto for i=1:Nza-1 if (za_ad(i)<=espesor_e)&&(espesor_e<=za_ad(i+1)) posicion_z=i; break end end %Condición de contorno horizontales sobre el electrodo visto % for i=1:length(ix_elect) % I=(ix_elect(i)-1)*Nza+posicion_z; % Aa(I,:)=Dza(I,:); ba(I,1)=0; % end % Condición de contorno verticales sobre el electrodo visto % for j=1:(posicion_z-1) % I=(ix_elect(1)-1)*Nza+j; %Pared izquierda % Aa(I,:)=Dx(I,:); ba(I,1)=0; % % I=(ix_elect(length(ix_elect))-1)*Nza+j; %Pared derecha % Aa(I,:)=Dx(I,:); ba(I,1)=0; % end % Condición de contorno sobre todo el electrodo for q=1:posicion_z for i=1:length(ix_elect) I=(ix_elect(i)-1)*Nza+q; Aa(I,:)=0; Aa(I,I)=1; ba(I,1)=0; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% DOMINIO DIELÉCTRICO %%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % MATRICES DERIVADAS [dzb,dzb2,~,~]=Derivs(zb_ad,Nzb); Dxb=sparse(kron(dx,speye(Nzb))) ; Dxb2=sparse(kron(dx2,speye(Nzb))); Dzb=sparse(kron(speye(Nx),dzb)); Dzb2=sparse(kron(speye(Nx),dzb2)); DLb=sparse(Dxb2+Dzb2) ; % PERMITIVIDAD RELATIVA DEL MEDIO eps_b=eps_kapton*ones(Ntb,1); for i=1:Nx eps_b((i-1)*Nzb+Nzb,1)=eps_mean; end % MATRIZ CON LA ECUACIÓN A RESOLVER Ab=sparse(eps_b.*DLb); % TÉRMINO INDEPENDIENTE DEL SISTEMA Ab·rhocb=bb bb(1:Ntb,1)=0; % CONDICIONES DE CONTORNO % 1) Contornos verticales: dphidxb=0 for j=2:Nzb-1 I=j; Ab(I,:)=0; Ab(I,I)=1; bb(I,1)=0; I=(Nx-1)*Nzb+j; Ab(I,:)=0; Ab(I,I)=1; bb(I,1)=0; end % 2) Contorno horizontal inferior; for i=2:(Nx-1) I=(i-1)*Nzb+1; Ab(I,:)=0; Ab(I,I)=1; bb(I,1)=0; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% ENTREFASE %%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % CONDICIONES DE CONTORNO Bb=zeros(Ntb,Nta); %Definición de la matriz Bb Ba=zeros(Nta,Ntb); %Definición de la matriz Ba j=1; for i=1:Nx Ia=(i-1)*Nza+1; Ib=(i-1)*Nzb+Nzb; if (x_ad(i)-ini_eo/L)*(fin_eo/L-x_ad(i))>=0 70 70 forcea=sqrt(forcexa.^2+forceza.^2); forceb=sqrt(forcexb.^2+forcezb.^2); % CREACIÓN MATRICES – Para graficar posteriormente for j=1:Nzb xmat(1:Nx,j)=x_ad(1:Nx)'; zmat(1:Nx,j)=zb_ad(j); phimat(1:Nx,j)=phib(((1:Nx)-1)*Nzb+j,1); forcexmat(1:Nx,j)=forcexb(((1:Nx)-1)*Nzb+j,1); forcezmat(1:Nx,j)=forcezb(((1:Nx)-1)*Nzb+j,1); forcemat(1:Nx,j)=forceb(((1:Nx)-1)*Nzb+j,1); Emat(1:Nx,j)=Eb(((1:Nx)-1)*Nzb+j,1); rhocmat(1:Nx,j)=rhocb(((1:Nx)-1)*Nzb+j,1); Exmat(1:Nx,j)=Exb(((1:Nx)-1)*Nzb+j,1); Ezmat(1:Nx,j)=Ezb(((1:Nx)-1)*Nzb+j,1); end for j=2:Nza xmat(1:Nx,j+Nzb-1)=x_ad(1:Nx)'; zmat(1:Nx,j+Nzb-1)=za_ad(j); phimat(1:Nx,j+Nzb-1)=phia(((1:Nx)-1)*Nza+j,1); forcexmat(1:Nx,j+Nzb-1)=forcexa(((1:Nx)-1)*Nza+j,1); forcezmat(1:Nx,j+Nzb-1)=forceza(((1:Nx)-1)*Nza+j,1); forcemat(1:Nx,j+Nzb-1)=forcea(((1:Nx)-1)*Nza+j,1); Emat(1:Nx,j+Nzb-1)=Ea(((1:Nx)-1)*Nza+j,1); rhocmat(1:Nx,j+Nzb-1)=rhoca(((1:Nx)-1)*Nza+j,1); Exmat(1:Nx,j+Nzb-1)=Exa(((1:Nx)-1)*Nza+j,1); Ezmat(1:Nx,j+Nzb-1)=Eza(((1:Nx)-1)*Nza+j,1); end % VARIABLES filename='Potencial.mat'; save(filename, 'rhoca',…); end En el caso del cálculo del potencial el proceso es análogo, se cargan los parámetros necesarios, se crean las matrices de derivadas, se define la matriz del sistema y el término independiente, se imponen las condiciones de contorno y finalmente se resuelve el sistema. De nuevo se tiene en cuenta el espesor del electrodo visto. Además, tras resolver el sistema se obtiene el potencial eléctrico y a continuación se calculan el campo eléctrico en cada punto, así como la fuerza electrohidrodinámica generada, término que posteriormente se incluye como motor del movimiento en el problema fluido. Problema fluido function Problema_fluido close all; clc; % PARÁMETROS [Nx,Nza,~,Nta,~,x_ad,za_ad,~,L]=parametros_mallado; [dtau,u0,Re,Dc]=parametros_fluido; % CARGA DATOS DE POTENCIAL carga=load ('Potencial.mat'); rhoc=carga.rhoca; Ex=carga.Exa; Ez=carga.Eza; Fx=carga.forcexa; Fz=carga.forceza; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% DOMINIO FLUIDO %%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Matrices derivadas [dxa, dxa2,dxa3,dxa4]=Derivs(x_ad,Nx); %X se malla con diferencias finitas [dza,dza2,dza3,dza4]=Derivs(za_ad,Nza); %Z se malla con diferencias finitas %[dza]=dCheby(za_ad,Nza); %Z se malla con Chebyshev Dxa =sparse(kron(dxa,speye(Nza))); Dxa2=sparse(kron(dxa2,speye(Nza))); Dxa3=sparse(kron(dxa3,speye(Nza))); Dxa4=sparse(kron(dxa4,speye(Nza))); Dza=sparse(kron(speye(Nx),dza)); Dza2=sparse(kron(speye(Nx),dza2)); Dza3=sparse(kron(speye(Nx),dza3)); Dza4=sparse(kron(speye(Nx),dza4)); DL=sparse(Dxa2+Dza2); %Matriz Laplaciano DL2=sparse(Dxa4+Dza4+2*Dxa2*Dza2); %Matriz Laplaciano^2 DxDL=sparse(Dxa3+Dxa*Dza2); DzDL=sparse(Dza*Dxa2+Dza3); % DEFINICIÓN DEL PROBLEMA % MATRIZ Lap-1/Re*Lap2*dt DEL SISTEMA A_psi*psi_n=b_psi APSI=sparse(DL-DL2/Re*dtau); % CONDICIONES DE CONTORNO % Factor de condiciones de contorno FBC_PSI=ones(Nta,1); % Condiciones que hacen Psi igual a cero en contornos horizontales: 72 72 APSI_tr=APSI'; for i=2:(Nx-1) % 1) Contorno inferior. Importante tener en cuenta que este contorno %inferior coincidirá con el contorno intermedio, donde se encuentra el %electrodo expuesto I=(i-1)*Nza+1; %psi=0 APSI_tr(:,I)=0; APSI_tr(I,I)=1; FBC_PSI(I,1)=0; % 2) Contorno superior I=(i-1)*Nza+Nza; %psi=0 APSI_tr(:,I)=0; APSI_tr(I,I)=1; FBC_PSI(I,1)=0; end % Condiciones que hacen psi igual a cero en contornos verticales: for j=1:Nza % 3) Contorno izquierdo I=j; %psi=0 APSI_tr(:,I)=0; APSI_tr(I,I)=1; FBC_PSI(I,1)=0; % 4) Contornos derecho I=(Nx-1)*Nza+j; %psi=0 APSI_tr(:,I)=0; APSI_tr(I,I)=1; FBC_PSI(I,1)=0; end %Falta por aplicar velocidades nulas en las paredes. Para ello, se usan las %filas correspondientes a los subcontornos de abajo. % (1-2) Contornos horizontales dpsi/dz=0 for i=3:(Nx-2) K=(i-1)*Nza+2; APSI_tr(:,K)=Dza(:,K-1); FBC_PSI(K,1)=0; K=(i-1)*Nza+Nza-1; APSI_tr(:,K)=Dza(:,K+1); FBC_PSI(K,1)=0; End % (3-4) Contornos verticales dpsi/dx=0 for j=2:(Nza-1) K=(2-1)*Nza+j; APSI_tr(:,K)=Dxa(:,K-Nza); FBC_PSI(K,1)=0; K=(Nx-2)*Nza+j; APSI_tr(:,K)=Dxa(:,K+Nza); FBC_PSI(K,1)=0; end APSI=APSI_tr'; % RESOLUCIÓN DEL SISTEMA: FUNCIÓN DE CORRIENTE % Descomposición LU [Lu,U,P]=lu(APSI); PpL=P'*Lu; for i=1:Nx for j=1:Nza I=(i-1)*Nza+j; Psinm1(I,1)=0; end end figure Ntiempo=200; for n=1:Ntiempo tau=n*dtau; ft=sin(2*pi*4500*L*tau/u0); disp(n) bLinPsi(1:Nta,1)=DL*Psinm1; %Término lineal bNLPsi(1:Nta,1)=(+(Dza*Psinm1).*(DxDL*Psinm1)- (Dxa*Psinm1).*(DzDL*Psinm1))*dtau; %Término no lineal bElectPsi=Dc*(Ex.*(Dza*rhoc)-Ez.*(Dxa*rhoc))*ft.^2*dtau; %Término del campo eléctrico bpsi=bLinPsi-bNLPsi+bElectPsi; bsyst(1:Nta,1)=FBC_PSI.*bpsi; aux=PpL\bsyst; Psi=U\aux; % CAMPO DE VELOCIDADES vx=Dza*Psi; vz=-Dxa*Psi; V=sqrt(vx.^2+vz.^2); 74 74 % CREACIÓN DE MATRICES for j=1:Nza xmat(1:Nx,j)=x_ad(1:Nx)'; zmat(1:Nx,j)=za_ad(j); psimat(1:Nx,j)=Psi(((1:Nx)-1)*Nza+j,1); vxmat(1:Nx,j)=vx(((1:Nx)-1)*Nza+j,1); vzmat(1:Nx,j)=vz(((1:Nx)-1)*Nza+j,1); Vmat(1:Nx,j)=V(((1:Nx)-1)*Nza+j,1); end % GRÁFICOS if n==5*round(n/5) %Graficamos las línes de corriente subplot(1,2,1), contour (xmat, zmat, psimat); %subplot(1,3,2),quiver(xmat,zmat,vxmat,vzmat); subplot(1,2,2),contourf(xmat,zmat,Vmat), axis([-2 10 -1 1.5]); pause(0.01) hold off end % ACTUALIZACIÓN Psinm1=Psi; end end Cálculo de matrices derivadas basadas en Diferencias Finitas function[df d2f d3f d4f]=Derivs(x, N) df=sparse(N,N); d2f=sparse(N,N); 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 Cálculo de las matrices derivadas según Chebyshev 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, 76 76 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 3 ANÁLISIS DE RESULTADOS El código numérico implementado a lo largo de este estudio tiene gran flexibilidad, pudiendo llevarse a cabo una gran serie de investigaciones con él. El primer tema a abordar será la validación del mismo, comparando los resultados arrojados por el código en el problema eléctrico con los resultados de la literatura, la validación del problema fluido se lleva a cabo en [9]. Tras comprobar la fiabilidad del código se analiza la influencia de distintos parámetros de diseño del actuador sobre el problema, tratando con especial atención su efecto sobre la fuerza electrohidrodinámica. Por otro lado, se implementa el problema con condiciones de contorno periódicas y se comprueba el efecto de estas sobre la densidad de carga, el potencial debido a campos externos y el campo eléctrico. Por último, se lleva a cabo una estimación del coste computacional, registrando el tiempo empleado por el código para resolver el problema. 3.1 Validación del código El código se divide en tres módulos, cada uno asociado a uno de los problemas fundamentales que han de resolverse: el cálculo de la densidad de carga, la obtención del potencial eléctrico debido a cargas externas y la resolución del problema fluido incluyendo la fuerza electrohidrodinámica en las ecuaciones de Navier-Stokes. A continuación, se cotejarán los resultados obtenidos en cada uno problemas con la literatura, con el fin de comprobar que los resultados que arroja el código implementado son veraces y de confianza. Parámetros Con el fin de poder comparar los resultados de forma objetiva, se eligen los mismos parámetros que los utilizados en el artículo [15].  Densidad de carga máxima, 𝜌𝑐,𝑚𝑎𝑥=0,0008𝐶 𝑚3  Amplitud AC, 𝜙𝑚𝑎𝑥=5,0 𝑘𝑉  Frecuencia AC, 𝑓=4,5 𝑘𝐻𝑧  Longitud de Debye, 𝜆𝐷=0,001 𝑚  Densidad del aire, 𝜌=1,225𝑘𝑔 𝑚3  Viscosidad dinámica del aire, 𝜇=1,85·10−5𝑘𝑔 𝑚·𝑠  Longitud electrodo visto =0,01 𝑚  Longitud electrodo embebido =0,01 𝑚  Longitud de la separación (gap) entre electrodos =0,0005 𝑚  Posición en Z del electrodo embebido =0,000125 𝑚  Espesor de los electrodos =0,0001 𝑚 Densidad de carga El funcionamiento del actuador de plasma se basa en su excitación mediante una corriente alterna de alto voltaje, que provoca que el aire en el entorno del dispositivo quede débilmente ionizado. Ese aire ionizado en presencia del gradiente de campo eléctrico entre los electrodos genera una fuerza volumétrica sobre el flujo, la cual induce componentes de velocidad sobre este. La densidad de carga mide la cantidad de carga concentrada en el dominio, que como se comprobará a continuación, se agrupa prácticamente en su totalidad sobre el electrodo embebido. 78 78 La densidad de carga obtenida en todo el dominio fluido aparece reflejada en la siguiente figura: Ilustración 37 – Densidad de carga obtenida con el código de Matlab. Código sencillo. Ilustración 38 - Densidad de carga obtenida con el código de Matlab. Código complejo. Ilustración 39 – Densidad de carga obtenida en [15] usando GHOST (Código desarrollado por la Universidad de Kentucky). Como se puede comprobar, las tres figuras presentan la misma escala de color y la misma distribución espacial. Los resultados son bastante parecidos, la distribución de densidad de carga presenta valores casi idénticos en todas las gráficas y una geometría similar. Cabe mencionar que Suzen et Al. en [15] utiliza un código tipo caja negra para resolver el problema, mientras que en nuestro caso se lleva a cabo con un método completamente casero. Por consiguiente, sus códigos presentan mayor complejidad y aproximan mejor el fenómeno físico. Se puede observar cómo el espesor de los electrodos presenta cierta influencia sobre la distribución de densidad de carga sobre el dieléctrico. En la literatura con la que se compara se ha considerado el problema de espesor, mientras que, en uno de los códigos desarrollados en este trabajo se ha supuesto despreciable por motivos de simplicidad y en otro no. A la vista está que las hipótesis simplificadoras tomadas modifican ligeramente el resultado, pero se obtienen datos que en primera aproximación son muy similares y requieren menor coste computacional, lo cual se traduce en tiempo. Potencial eléctrico debido a cargas externas Se recuerda que en el caso de gases débilmente ionizados el potencial eléctrico se puede dividir en dos términos, Φ=𝜙+𝜑 El primero de ellos es debido al campo eléctrico externo que se aplica sobre los electrodos y no tiene un impacto considerable en la distribución de densidad de carga. El segundo es el potencial debido a la densidad de carga neta del plasma y es el que determina cómo es la distribución de esta a lo largo del dominio. A continuación, se comparará el potencial debido a cargas externas, 𝜙, con los resultados obtenidos por Suzen et Al [15] y por Seth en [20]. Esta magnitud permite el posterior cálculo del campo eléctrico y junto con la densidad de carga de la fuerza electrohidrodinámica. Ilustración 40 – Potencial eléctrico obtenido usando Matlab. Código sencillo. 86 86 Ilustración 52 – Región para el cálculo de la fuerza media. En la Ilustración 53 se representa la fuerza máxima total, junto con la fuerza máxima en la dirección x y z. Se puede observar cómo la componente longitudinal de la fuerza es dominante sobre la vertical, lo cual es un hecho típico en los actuadores DBD. A medida que el espesor del dieléctrico disminuye los electrodos se encuentran más cerca y el campo eléctrico se hace más intenso y, por consiguiente, la fuerza total máxima se ve incrementada. Ilustración 53 – Fuerza máxima en función del espesor de dieléctrico. Además, la fuerza longitudinal es muy cercana a la fuerza total, de ahí que el actuador de plasma inyecte una mayor componente de velocidad en la dirección x que en z. Por otro lado, puede parecer que el hecho de variar la altura de la capa dieléctrica no presenta gran repercusión sobre la magnitud de la fuerza máxima inducida en la componente z, que aparece al final del electrodo visto y que mantiene en un valor prácticamente constante en la figura anterior, esto se debe a que la escala de esta solicitación es mucho menor. En la Ilustración 54 puede observarse cómo al variar la altura de la capa dieléctrica desde 5 mm hasta 0,127 mm la magnitud de la fuerza vertical se triplica. Ilustración 54 – Fuerza máxima en z en función del espesor de dieléctrico. Como se mencionó anteriormente, la disminución del espesor del dieléctrico implica que el campo eléctrico es más intenso y con ello que la fuerza también lo es. Sin embargo, los resultados anteriores no permiten obtener una idea concisa sobre cómo es la distribución de las fuerzas en el entorno de los electrodos y cúal es el área de influencia de esta. Para ello, se estudia la fuerza longitudinal media dentro de la región fijada en la Ilustración 52 a medida que el espesor de la capa aislante varia. Ilustración 55 – Fuerza media en x en función del espesor del dieléctrico. Como se puede apreciar, la fuerza media máxima ocurre cuando el espesor del dieléctrico está en torno a 1 mm. Para espesores menores, aunque la fuerza máxima es cada vez mayor, en general su zona de influencia disminuye. Mientras que, para espesores mayores, tanto la fuerza máxima como el área donde la solicitación se concentra decrecen con la altura del dieléctrico. Relación Fuerza – Gap entre electrodos El gap es la distancia longitudinal entre el extremo derecho del electrodo visto y el izquierdo del electrodo embebido, un parámetro geométrico a tener en cuenta al diseñar el actuador puesto que, de nuevo, la distancia entre electrodos se ve modificada y con ella el campo eléctrico y la fuerza generada. 88 88 En el análisis anterior se estableció un gap fijo de 0,5 mm y se estudió el efecto de aumentar la altura del dieléctrico, obteniendo un espesor óptimo de 1 mm que maximizaba el área de actuación de la fuerza electrohidrodinámica. Siguiendo el mismo procedimiento se fija el espesor de dieléctrico al óptimo y se estudia el comportamiento del actuador ante distintos valores del gap, que varían entre 0 y 5 mm. Nuevamente, se evalua la evolución de la fuerza máxima y la fuerza media para analizar el efecto del gap sobre el actuador. La fuerza máxima viene representada en la Ilustración 56, se puede observar cómo esta aumenta a medida que el gap tiende a 0. Por tanto, de nuevo se comprueba cómo al disminuir la distancia entre electrodos la fuerza máxima se ve incrementada. Ilustración 56 – Fuerza máxima en función del gap. Además en la Ilustración 57, para el caso de la fuerza media en la dirección x, que es la componente predominante en el campo de fuerzas, se verifica también que aumenta a medida que la distancia entre los electrodos tiende a cero. Ilustración 57 – Fuerza media en x en función del gap. Relación Fuerza – Voltaje aplicado La relación entre la fuerza eléctrica producida por el actuador y el voltaje aplicado sobre los electrodos es lineal. La fuerza calculada con el código elaborado está adimensionalizada con el propio voltaje impuesto sobre el actuador, por tanto, ambas magnitudes son directamente proporcionales. Ilustración 58 – Relación fuerza máxima – voltaje. De manera análoga, esa relación lineal se transmite a la velocidad inducida sobre el campo fluido, como se puede comprobar en la publicación de Dugo, P [9], de la que se toma la Ilustración 59. Ilustración 59 – Relación velocidad inducida – voltaje [9]. 3.3 Influencia de los parámetros experimentales del modelo de SH sobre el comportamiento del actuador El modelo SH es un modelo macroscópico, que a diferencia de los modelos microscópicos que calculan la trayectoria de cada una de las partículas cargadas que se acumulan sobre la superficie del electrodo embebido, considera a la densidad de carga como un todo, un ente único y que se impone en función de resultados experimentales. Uno de los dos parámetros que definen la fuerza generada por el actuador es la propia distribución de densidad de carga y el otro es el campo eléctrico, como se puede comprobar en la expresión (12). El hecho de que la densidad de carga se calibre para este modelo de forma experimental implica que dicha distribución es sólo válida para el caso establecido durante el experimento, y necesitaría ser calibrada para las distintas configuraciones que se impongan sobre el actuador. 90 90 Por otro lado, el modelo SH también impone el valor de la longitud de Debye, parámetro que de nuevo se obtiene mediante experimentación y que influye directamente sobre cómo se concentran las cargas sobre la pared del electrodo enterrado, afectando así a la propia fuerza eléctrica y al campo de velocidades. Así, tanto la longitud de Debye como la densidad de carga máxima son parámetros a calibrar para cada posible configuración del actuador, por ello resulta interesante estudiar la evolución del problema a medida que ambos parámetros cambian. Cabe mencionar que el potencial asociado a campos externos, y por consiguiente el campo eléctrico, no se ve alterado por dichas magnitudes al no estar presentes en la ecuación que gobierna su comportamiento (15). Influencia de la longitud de Debye sobre el problema eléctrico La longitud de Debye es una magnitud característica del plasma, que mide el alcance del efecto electrostático de una partícula cargada dentro del propio plasma. El alcance es limitado porque los electrones tienden a agruparse en torno a las partículas cargadas positivamente y acaban apantallando el efecto electrostático de estas, que sólo será detectado por los electrones que la rodean. Por consiguiente, el campo eléctrico creado por un ion positivo está confinado dentro de una región que se extiende a cierta distancia alrededor de la partícula, dicha distancia es la longitud de Debye. En la siguiente imagen se representa el significado físico de este parámetro. Ilustración 60 – Fenómeno de apantallamiento. Tanto la longitud de Debye como la densidad de carga máxima del plasma son parámetros desconocidos en el modelo SH, son magnitudes que cambian según las condiciones físicas del plasma, sobre todo con la temperatura, y que son necesario calibrar con experimentación. Anteriormente en esta investigación se han tomado los valores de estos parámetros basados en medidas experimentales impuestos en el artículo de Suzen et Al. [15]. Sin embargo, el valor teórico de la longitud de Debye para fenómenos con descarga de gases a temperatura atmosférica es del orden de 10−4 m, como se puede encontrar en Blandford et al. (2004) [22] o en Roth (1995) [23]. Se pueden encontrar numerosos artículos donde se estudia el efecto de variar 𝜆𝐷, por ejemplo la publicación de Ibrahim et Al. Donde se cuantifica el resultado de estas variaciones sobre la distribución de densidad de carga y sobre la velocidad inducida en el flujo [21]. En dicha investigación se comprueba que a medida que incrementa 𝜆𝐷 la componente X de las velocidades aumenta pero a la vez se pierde precisión en la velocidad en la componente Z. A continuación, se llevará a cabo un pequeño estudio con el modelo desarrollado en Matlab, donde se comprobará el efecto de variar la longitud de Debye entre 0.1 mm y 5 mm y cómo afecta esto a la distribución de densidad de carga y a la fuerza electrohidrodinámica. Este estudio se basa en el llevado a cabo por Ibrahim, I.H. en [21]. En la publicación [21] se desea investiga el efecto de variar la longitud de Debye sobre la densidad de carga y la fuerza justo en la proyección sobre el dieléctrico del extremo izquierdo del electrodo embebido, es decir, el punto donde la fuerza es máxima según las condiciones de contorno impuestas. Sin embargo, se considera que realizar medidas en dicho punto carece de sentido físico; en primer lugar porque el modelo numérico impone condiciones de contorno fijas, haciendo que la densidad de carga sobre ese punto sea siempre máxima e igual a 1, independientemente del valor de 𝜆𝐷 que se imponga, hecho que indica que los resultados presentados por Ibrahim, I.H. no son correctos. Por otro lado, estudiar la fuerza sobre la pared no tiene ningún tipo de interés físico puesto que en el problema fluido se imponen condiciones de no deslizamiento, haciendo que la velocidad del flujo en los puntos de la superficie del dieléctrico sea nula. De acuerdo con lo anteriormente expuesto, en este estudio se recoge la influencia de la longitud de Debye, tomando medidas a cierta distancia sobre la pared del dieléctrico, concretamente a 3·10−5 m. Esto permite cuantificar cómo varía la distribución espacial de la densidad de carga y también la evolución de la fuerza sobre un punto que está a una distancia normal a la pared fija y distinta de cero. Ilustración 61 – Densidad de carga adimensional en función de la longitud de Debye. En la imagen anterior se representa la densidad de carga adimensional sobre el extremo izquierdo del electrodo oculto y a 3·10−5 𝑚 sobre la alores pequeños de la longitud de Debye las cargas se acumulan muy cerca de la pared. superficie del dieléctrico. Se puede observar cómo la longitud de Debye influye en la propagación de la carga sobre el electrodo. Para 𝜆𝐷 del orden de 10−5 𝑚 la carga acumulada sobre el electrodo presenta una extensión insignificante, quedando concentrada casi en su totalidad por debajo del punto en el que estamos midiento. Por el contrario, a partir de longitudes del orden de 10−4 𝑚 la densidad de carga aumenta considerablemente, pasando de estar en 0.15 cuando la longitud de Debye era 10−5𝑚 a 0.77 cuando la longitud de Debye es 10−4 𝑚, esto indica que a medida que se ha aumentado 𝜆𝐷 el apantallamiento de las cargas necesita mayor distancia y la acumulación de cargas ocupa más espacio. A partir de ahí, la propagación espacial de la densidad de carga no es tan acusada, ya que el punto sobre el que se mide está contenido sobre la región donde se concentran las cargas. Por tanto, se puede concluir que para v En la siguiente imagen se representa la densidad de carga sobre la proyección a 3·10−5𝑚 sobre la pared dieléctrica del extremo izquierdo del electrodo embebido frente a distintos valores de la longitud de Debye, lo que permitirá comprender cómo es la relación entre ambas magnitudes. 92 92 Ilustración 62 – Densidad de carga máxima adimensional en función de la longitud de Debye. Efectivamente, se puede comprobar cómo para valores pequeños de la longitud de Debye la variación de la densidad de carga es muy brusca, mientras que para valores en torno a 10−4 m aparece una relación aproximadamente lineal entre el crecimiento de la densidad de carga y de la longitud de Debye. Según se explica en la publicación de Ibrahim, I.H. y Skote, M. [21], el notable crecimiento en la densidad de carga para longitudes de Debye pequeñas se debe a que en esta zona el valor de la densidad carga está condicionado a la permitividad relativa del dieléctrico, en nuestro caso Kapton, con 𝜀=2.7. Para valores de la longitud de Debye mayores, la densidad de carga pasa a estár condicionada por la permitividad del aire, que es 𝜀=1, haciendo que el crecimiento sea lineal con la longitud de Debye. Sin embargo, en este trabajo se rebate dicha explicación, el crecimiento en la densidad de carga se debe a que para longitudes de Debye menores de la altura sobre la que se está midiendo, 3·10−5𝑚, la carga se concentra por debajo del punto de medida y, por tanto, el valor de densidad de carga medido es pequeño. Una vez que 𝜆𝐷 aumenta y sobrepasa los 3·10−5 𝑚 se puede comprobar cómo la densidad de carga comienza a aumentar de forma lineal. Por último, se adjuntan una serie de imágenes de la densidad de carga para distintas longitudes de Debye. Se puede comprobar cómo para 𝜆𝐷=0.001, la carga se concentra hasta una altura de en torno a 2.5 mm. 𝜆𝐷=1·10−5𝑚 𝜆𝐷=5·10−5𝑚 𝜆𝐷=1·10−4𝑚 𝜆𝐷=2·10−4𝑚 𝜆𝐷=3·10−4𝑚 𝜆𝐷=4·10−4𝑚 𝜆𝐷=5·10−4𝑚 𝜆𝐷=6·10−4𝑚 𝜆𝐷=7·10−4𝑚 𝜆𝐷=8·10−4𝑚 𝜆𝐷=9·10−4𝑚 𝜆𝐷=1·10−3𝑚 La densidad de carga influye de manera directa sobre la fuerza generada por el actuador, así, repitiendo el análisis llevado a cabo anteriormente, pero esta vez con el fin de observar los efectos de variar la longitud de Debye sobre la fuerza se obtienen los siguientes resultados. 94 94 Ilustración 63 – Fuerza adimensional en función de la longitud de Debye. De manera análoga al caso de la densidad de carga, la fuerza presenta una variación brusca para longitudes pequeñas de Debye, lo cual es obvio, puesto que para 𝜆𝐷 ~ 10−5 m la carga se acumula casi en su totalidad por debajo del punto sobre el que tomamos la fuerza. A medida que aumenta la longitud de Debye su propagación espacial aumenta y la fuerza en dicha altura fija también lo hace. Una vez alcanzada una longitud del orden de 10−4 𝑚 la relación entre estas dos magnitudes comienza a ser aproximadamente lineal, muestra de ello es la siguiente figura que representa la fuerza adimensional, situada a 3·10−5𝑚 sobre la proyección en el dieléctrico del extremo izquierdo del electrodo embebido, en función de la longitud de Debye. Ilustración 64 – Fuerza adimensional máxima en función de la longitud de Debye. Al igual que con la densidad de carga se adjunta una serie de representaciones escalares del módulo de la fuerza para distintos valores de la longitud de Debye, que permite observar la relación entre ambas magnitudes. 𝜆𝐷=1·10−5𝑚 𝜆𝐷=5·10−5𝑚 𝜆𝐷=1·10−4𝑚 𝜆𝐷=2·10−4𝑚 𝜆𝐷=3·10−4𝑚 𝜆𝐷=4·10−4𝑚 𝜆𝐷=5·10−4𝑚 𝜆𝐷=6·10−4𝑚 𝜆𝐷=7·10−4𝑚 𝜆𝐷=8·10−4𝑚 𝜆𝐷=9·10−4𝑚 𝜆𝐷=1·10−3𝑚 Por último, el efecto sobre el problema fluido al aumentar la longitud de Debye es el aumento de sustancial de la velocidad inducida, además el punto de máxima velocidad inducida se desplaza aguas abajo. En la Ilustración 65 se plasma lo anteriormente comentado. El estudio detallado de la influencia de la longitud de 102 102 se ve modificado, la densidad de carga sigue acumulándose de la misma forma sobre cada uno de los electrodos embebidos. Por tanto, la fuerza eléctrohidrodinámica lejos del comienzo de la proyección vertical electrodo embebido es nula y hace que los actuadores se comporten de forma independiente. No obstante, cabe recordar que la distribución de densidad de carga es un parámetro experimental impuesto en el modelo SH y que no contempla ni la variación temporal ni con el voltaje de esta. Por tanto, aunque con la distribución semi – Gaussiana empleada los actuadores se comporten de forma independiente, sería necesario realizar un estudio con una distribución más precisa y que modele de forma más realista el proceso de acumulación de las cargas sobre el electrodo oculto. 5 LÍNEAS FUTURAS Los actuadores de plasma son una tecnología en vigente desarrollo, debido a esto, los posibles campos de investigación sobre los que se puede trabajar son innumerables. A lo largo del desarrollo del trabajo se han encontrado ciertos puntos que podrían ser tratados con mayor detenimiento, a continuación, se recoge una serie de puntos cuyo estudio resultaría de interés de cara a proyectos futuros:  Estudio de la distribución de la densidad de carga sobre el electrodo embebido, el modelo de SH establece una distribución semi Gaussiana basada en medidas experimentales de la luz emitida por el plasma. Sin embargo, no considera la variación temporal ni con el voltaje de la densidad de carga en el actuador. Esto podría tener consecuencias en la relación entre actuadores adyacentes al imponer condiciones de contorno periódicas.  Ampliar la flexibilidad del método numérico elaborado en este proyecto adaptándolo a dominios con distinta geometría, así se podría cuantificar el efecto del actuador sobre superficies curvadas, como las de un perfil aerodinámico.  Estudiar distintas alternativas a las condiciones de contorno impuestas sobre la entrefase entre el medio fluido y el medio dieléctrico. En la literatura viene poca información sobre por qué se imponen cada una de las condiciones, no existe un concenso que establezca que esas son las condiciones de contorno adecuadas, cada estudio utiliza aquellas que considere que modelan mejor la física del problema.  Construcción de un actuador de plasma, debido a su arquitectura simple es un proyecto sencillo y barato, que permitiría el posterior estudio experimental de estos dispositivos. Sí requiere más trabajo la preparación del montaje eléctrico asociado, es decir, un generador de corriente de alto voltaje que permita conseguir los requisitos del actuador y que permita variar los parámetros eléctricos.  Estudio experimental en túnel de viento de un actuador de plasma sencillo, comprobando su eficacia a la hora de reenergizar la corriente sobre una superficie fluida. Además, se podría llevar a cabo una validación experimental del código desarrollado en este trabajo, midiendo la fuerza electrohidrodinámica y la velocidad que es capaz de inducir el actuador sobre el flujo.  Desarrollo de distintas estrategias de control para actuadores DBD que permita su utilización como dispositivo para reducir la resistencia aerodinámica, reenergizar la capa límite, mantener la corriente adherida o retrasar su desprendimiento, sustituir a los mandos de control de vuelo móviles, o ayudar en la refrigeración de ciertas superficies. 104 104 BIBLIOGRAFÍA [1]. Roth JR, Sherman DM, Wilkinson SP. (1998). Boundary layer flow control with a one atmosphere uniform glow discharge surface plasma. AIAA Meeting (Reno, USA, January 1998) paper #98-0328. [2]. Roth JR. (1998). Electrohydrodynamically induced airflow in a one atmosphere uniform glow discharge surface plasma. 25th IEEE Int. Conf. Plasma Science (Raleigh, USA). [3]. Benard N, Jolibois J, Moreau E. (2009). Lift and drag performances of an axisymmetric airfoil controlled by plasma actuator. J Electrostat 67:133–139. [4]. Corke TC, Post ML, Orlov DM. (2007). SDBD plasma enhanced aerodynamics: concepts, optimization and applications. J Progr Aerosp Sci 43:193–217. [5]. Sung Y, Kim W, Mungal MG, Cappelli MA. (2006). Aerodynamic modification of flow over bluff objects by plasma actuation. Exp Fluids 41:479–486. [6]. Thomas FO, Kozlov A, Corke TC. (2006). Plasma actuators for bluff body flow control. AIAA Paper 2006-2845. [7]. Li Y, Zhang X, Huang X. (2010). The use of plasma actuators for bluff body broadband noise control. Exp Fluids 49:367–377. [8]. Grundmann S, Tropea C. (2008). Active cancellation of artificially introduced Tollmien–Schlichting waves using plasma actuators. Exp Fluids 44:795–806. [9]. Dugo, P. (2019). Modelado y estudio fluidomecánico de los actuadores de plasma SDBD. Trabajo de Fin de Grado. Universidad de Sevilla. Sevilla. [10]. Corke, T., Jumper, E., Post, M., Orlov, D., McLaughlin, T. (2002). Application of Weakly Ionized Plasmas as Wing Flow Control Devices. AIM-2002-0350. [11]. Hultgren, L.S., Ashpis, D.E., 2003, Demonstration of Separation Delay with Glow-Discharge Plasma Actuators. AIAA-2003-1025. Reno, NV, Jan. 2003. [12]. Massines F, Rabehi A, Decomps P, Gadri RB, SeÌgur P and Mayoux C (1998). Experimental and theoretical study of a glow discharge at atmospheric pressure controlled by dielectric barrier. Journal of Applied Physics 83(6): 2950-2957. [13]. Boeuf JP and Pitchford LC (2005). Electrohydrodynamic force and aerodynamic flow acceleration in surface dielectric barrier discharge. Journal of Applied Physics 97(10): 1-10. [14]. Rouffet J-B (2008). Plasma Actuator Influence on Air Flow. COMSOL Conference 2008, Hannover, COMSOL. [15]. Suzen, Y.B., Huang, P.G., Jacob, J.D., Ashpis, D.E. (2005). Numerical Simulations of Plasma Based 106 106 Flow Control Applications. AIAA 2005-4633. [16]. Henry DC (1931). The Cataphoresis of Suspended Particles. Part I. The Equation of Cataphoresis. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 133(821): 106 - 129. [17]. C. L. Enloe, Thomas E. McLaughlin, Robert D. VanDyken, and John C. Fischer. (2004). Mechanisms and responses of a single dielectric barrier plasma actuator: Geometric effects. 42nd AIAA Aerospace Sciences Meeting and Exhibit. [18]. Ruiz Contreras, P.J. (2016). El método de colocación para el problema de convección de Rayleigh- Bénard rectangulares (Proyecto Mecánica de Fluidos Avanzada). Trabajo de Fin de Grado. Universidad de Sevilla, Sevilla. [19]. U.K. Seth, P. Traoré, F.J. Duran-Olivencia, E. Moreau, P.A. Vazquez. (2018). Parametric study of a DBD plasma actuation based on the Suzen-Huang model. Journal of Electrostatics. 93 (pp. 1-9). [20]. U.K. Seth, P. Traoré, F.J. Duran-Olivencia, E. Moreau, P.A. Vazquez. (2018). Application to the numerical simulation of various typical Electro-hydro-dynamic flows. Tesis doctoral. Universidad de Poitiers, Poitiers. [21]. I. H. Ibrahim and M. Skote. (2013). “Effects of the scalar parameters in the Suzen-Huang model on plasma actuator characteristics,” International Journal of Numerical Methods for Heat and Fluid Flow, August 2013. [22]. Blandford RD and Thorne K (2004). Applications of Classical Physics. California, USA, Caltech: p. 980. [23]. Roth JR (1995). DC Electrical Glow Discharge in Gases. Industrial Plasma Engineering Volume 1: Principles, Institute of Physics Publishing 1: Principles: 284-351. [24]. Bouchmal. A. (2011). Modelling of Dielectric-Barrier Discharge Actuator. Thesis. Delft University of Technology. [25]. Orlov, D.M. Corke, T.C. (2006). Electric circuit model for aerodynamic plasma actuator. AIAA Aerospace Sciences Meeting and Exhibit, 44. [26]. Chaves, W. Garcia, F. (2008). Espacio Virtual de la Física, Plasma. Volumen 6, Número 1. Instituto Tecnológico de Costa Rica. [27]. J. R. Reitz, F. J. Milford y R. W. Christy. (1996). Fundamentos de la teoría electromagnética, 4º ed, Cap. 14. Addison-Wesley Iberoamericana. [28]. www.temasek-labs.nus.edu.sg [29]. www.pbs.org/wgbh/nova/article/plasma-air-control/ [30]. A S Taleghani, A Shadaram, M Mirzaei, S Abdolahipour. (2018). Parametric study of a plasma actuator at unsteady actuation by measurements of the induced flow velocity for flow control. Journal of the Brazilian Society of Mechanical Sciences and Engineering. 40 (4) (pp. 1-13). [31]. Shieh, T.H. (2017). Study of influencing characteristics on boundary-layer separation controlled by using DBD plasma actuator with modified model. International Journal of Heat and Mass Transfer. 113 (pp.1212- 1233). [32]. Zhang, X., Li, H., Choi, K.S., Song, L. (2017). Coherent structures induced by dielectric barrier discharge plasma actuator. Modern Physics Letters B. 31 (32). [33]. Erfani, R., Zare-Behtash, H., Hale, C., Kontis, K. (2015). Development of DBD plasma actuators: The double encapsulated electrode. Acta Astronautica. 109 (pp. 132-143). [34]. Balsera Barquero, Mª. S. (2018). Design of a Plasma Actuator Based on Dielectric Barrier Discharge (DBD) (Master’s Thesis). Universidad de Sevilla, Sevilla. [35]. Pérez Formigó, M. (2010). Capítulo 3: Descargas eléctricas en gases: concepto, clasificación y fenómenos asociados (pp. 12-26). Modelado de descargas eléctricas en gases mediante redes neuronales. Trabajo Fin de Master. Universidad de Sevilla, Sevilla. [36]. Babou, Y., Nieto Martin, E., Fajardo Peña, P. (2017). Simple body force model for Dielectric Barrier Discharge plasma actuator. European Conference for Aeronautics and Aerospace Sciencies (EUCASS). [37]. Méndez Montoro de Damas, A. (2018). Plasmas de descarga de barrera dieléctrica con empaquetamiento ferroeléctrico a presión atmosférica aplicados a la producción de hidrógeno, amoníaco y a la descontaminación del aire. (Tesis Doctoral Inédita). Universidad de Sevilla, Sevilla. [38]. Morales Pérez, C. (2016). 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. [39]. Imran Halimi, I., Skote, M. (2011). Boundary condition modifications of the Suzen-Huang plasma actuator model. International Journal of Flow Control. [40]. Vázquez Roy, J.L. (2009). Apuntes de Campos Electromagnéticos Curso 08-09. Universidad Carlos III de Madrid, Madrid. [41]. T. Brauner, S. Laizet, N. Benard and E. Moreau. (2016). “Modelling of Dielectric Barrier Discharge Plasma Acutuators for Direct Numerical Simulatons,” AIAA Flow control conference. June 2016. [42]. Pardo San Gil, R. Métodos numéricos en I.Q. Universidad Complutense de Madrid. Madrid. [43]. Miguel Pérez-Saborid, (2019). Notas sobre el cálculo general y eficiente de matrices de derivadas mediante los polinomios de Lagrange y sus aplicaciones. ETSI. Sin publicar. 108 108 GLOSARIO DBD: Dielectric Barrier Discharge. SDBD: Single Dielectric Barrier Discharge. AC: Corriente Alterna. Modelo SH: Modelo numérico Suzen-Huang. TNL: Término No Lineal.