Cálculo de capas límites en régimen hipersónico

Equation Chapter 1 Section 1 Proyecto Fin de Grado Grado en Ingeniería Aeroespacial Cálculo de Capas Límites en Régimen Hipersónico Autor: Ángel Ruy-Díaz Rojas Tutor: Miguel Pérez-Saborid Sánchez Pastor Dep. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Sevilla, 2017 iii Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Estudio de la Capas Límite en Régimen Hipersónico Autor: Ángel Ruy-Díaz Rojas Tutor: Miguel Pérez-Saborid Sánchez-Pastor Profesor titular Dep. de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2017 v Trabajo Fin de Grado: Estudio de la Capas Límite en Régimen Hipersónico Autor: Angel Ruy-Díaz Rojas 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, 2017 El Secretario del Tribunal vii A mi familia A mi tutor A Cristina ix Agradecimientos Este trabajo ha supuesto un gran esfuerzo, no solo para mi, sino para todos los que han estado a mí alrededor en esto meses que he dedicado. En primer lugar, me gustaría agradecer a mi familia el apoyo moral que se me ha dado, algo clave en los momentos más difíciles, y que ha permitido que pudiede seguir adelante. Por otra parte, alguien fundamental en el desarrollo del trabajo ha sido Miguel Pérez-Saborid, tutor de este trabajo. Sin su ayuda y las horas dedicadas por su parte no hubiera sido posible la realización y finalización de este trabajo. Además, me gustaría agradacer a mis compañeros no solo el apoyo moral sino el logístico, siendo generosos en ayudarme a ahorrar tiempo haciendo funcionar los progamas de MATLAB en sus ordenadores puediendo así ahorrar muchas horas de cálculo. Por último, pero no menos importante quería a agradecer a Cristina todo el apoyo que me ha ofrecido ya que ha estado donde debía estar y cuando debía para que no bajara los brazos y me rindiera. De corazón, muchas gracias a todos los mencionados y a los que habré olvidado mencionar porque sin el apoyo y la ayuda de todos y cada uno de vosotros, estre trabajo no hubiera sido posible. Angel Ruy-Díaz Rojas Sevilla, 2017 Figura 2-21 Esquema de perfil parabólico. 38 Figura 2-22 Coeficiente de presiones para perfil parabólico con 𝜎=0.2. 39 Figura 2-23 Coeficiente de presiones para perfil parabólico con σ=0.7 40 Figura 2-24 Coeficiente de presiones sobre distintos perfiles parabólicos a Mach 20 40 Figura 2-25 Coeficiente de resistencia de presión frente al número de Mach para distintos valores de 𝜎 41 Figura 2-26 Coeficiente de presiones sobre cilindro (M=20). 42 Figura 2-27 Coeficiente de resistencia en cilindro frente a Mach. 43 Figura 2-28 Placa plana con borde redondeado. 43 Figura 2-29 Coeficiente de presiones sobre placa plana con borde redondeado. 44 Figura 2-30 Coeficiente de resistencia de presión sobre placa plana con borde redondeado 44 Figura 2-31 Esquema de esfera. 45 Figura 2-32 Coeficiente de resistencia experimental de un cono y una esfera de [8]. 46 Figura 2-33 Coeficiente de resistencia aproximado para esfera y cono 47 Figura 3-1 Perfil de velocidad en capa límite. 50 Figura 3-2 Mallado uniforme con L puntos en el eje x y M puntos en la dirección y. 55 Figura 3-3 Mallado variable en x. 55 Figura 3-4 Perfil de velocidades en placa plana fria. 59 Figura 3-5 Perfiles de temperatura para placa plana con temperatura impuesta. 60 Figura 3-6 Perfil de velocidades en placa plana adiabática. 61 Figura 3-7 Perfil de temperaturas en placa plana adiabática. 62 Figura 3-8 Coeficiente de resistencia de fricción sobre placa plana laminar de [9] 63 Figura 3-9 Coeficiente de resistencia de fricción en una placa plana para diferentes relaciones de temperaturas. 64 Figura 3-10 Coeficiente de resistencia de fricción en perfiles elípticos 65 Figura 3-11 Evolución de perfil de velocidades a Mach 5 (izqda) y 25 (drcha). 65 Figura 3-12 Evolución de perfil de temperaturas a Mach 5 (izqda) y 25 (drcha). 66 Figura 3-13 Coeficiente de resistencia de fricción en perfiles parabólicos. 67 Figura 3-14 Evolución de perfil de velocidades a Mach 5 (izqda) y 25 (drcha) 67 Figura 3-15 Evolución de perfil de temperaturas a Mach 5 (izqda) y 25 (drcha). 68 Figura 3-16 Coeficiente de resistencia de fricción en perfiles elípticos de pared adiabática. 68 Figura 3-17 Coeficiente de resistencia de fricción para perfiles parabólicos de pared adiabática. 69 Figura 3-18 Distribución de transferencia de calor en cilndro (a la izquierda resultados numéricos, a la derecha resultados de la literatura). 71 Figura 3-19 Distribución de transferencia de calor en una placa plana con borde redondeado (a la izqda. resultados numéricos, a la drcha. resultados de la literatura). 72 Figura 3-20 Distribución de transferencia de calor para varios números de Mach a lo largo de un perfil elíptico con 𝜎=0.2 73 Figura 3-21 Distribución de transferencia de calor para diferentes perfiles elípticos a Mach 20. 74 Figura 3-22 Transferencia de calor en el punto de remanso. 74 Figura 3-23 Transferencia de calor total para perfiles elípticos a diferentes números de Mach y relaciones de aspecto. 75 xvii Figura 3-24 Variación del número de Nusselt en perfiles elípticos con 𝜎 y Mach. 76 Figura 3-25 Distribución de calor en la superficie para varios números de Mach. 76 Figura 3-26 Distribución de trasnferencia de calor para varias relaciones de aspecto. 77 Figura 3-27 Transferencia de calor en las inmediaciones del punto de remanso para perfiles elípticos 77 Figura 3-28 Transferencia total de calor en perfiles parabólicos 78 Figura 3-29 Dependencia del número de Nusselt con los diferentes parámetros. 79 Figura 3-30 Temperatura en perfil elíptico (𝜎=0.6) adiabático. 80 Figura 3-31 Temperatura en perfiles elípticos (M=21) adiabático. 80 Figura 3-32 Temperatura en la pared en perfile parabólico (𝜎=0.6) para diferentes números de Mach 81 Figura 3-33 Temperatura en la pared de difernetes perfiles parabólicos (𝑀∞=21). 81 1 1 INTRODUCCIÓN n el presente trabajo se va a hacer un estudio preliminar de la capa límite en régimen hipersónico. Sin embargo, para llegar hasta este punto es necesario responder antes a ciertas cuestiones: ¿Qué es el régimen hipersónico? ¿Es realmente relevante su estudio comoo para separarlo del régimen supersónico? 1.1. Régimen hipersónico De forma general, cuando se habla del régimen hipersónico, suele describirse como el régimen que se alcanza cuando el número de Mach es mayor que 5. También suelen mencionarse efectos de alta temperatura, reacciones químicas en el aire, etc. Por el contrario, no existe un salto cualitativo con el régimen supersónico como si ocurre cuando el número de Mach supera el valor unidad Sin embargo, una mejor definición de este régimen puede encontrarse en [1]. De aquí es de donde se obtiene la definición de régimen hipersónico que se seguirá en el presente trabajo. En él, se define como un régimen dentro del supersónico donde, al aumentar el número de Mach, se tornan relevantes varios de los siguientes efectos físicos: fina capa de aire tras la onda de choque, capa entrópica, interacción viscosa, flujo a alta temperatura y flujo de baja densidad. A continuación, se desarrollarán brevemente los distintos efectos. Volviendo a la cuestión del número de Mach al que el flujo puede considerarse hipersónico, en realidad no existe un número exacto. Se toma Mach igual a 5 como número intermedio de referencia porque es necesario hacer una distinción ya que los métodos para resolver las ecuaciones no son los mismos. Varios de los efectos son importantes a Mach menores que 5 y otros a velocidades mayores, por tanto, como primera aproximación, Mach 5 es un buen número, aunque 4 o 6 hubieran sido igualmente válidos. E 2 Estudio de la Capas Límite en Régimen Hipersónico Los diferentes efectos que suelen darse en el régimen hipersónico son los siguientes: -Fina capa de aire tras la onda de choque: Con el aumento del Mach de vuelo, la densidad tras las ondas de choque aumenta, por tanto, debido a la conservación de masa, la onda se mantiene muy cerca del cuerpo. Esto puede llegar a ser un problema a bajos números de Reynolds donde la capa viscosa puede aumentar mucho su grosor. Por el contrario, es de gran ayuda a valores altos de este número, donde al ser la capa prácticamente no viscosa permite usar la aproximación llamada thin shock layer theory [1]. Esta teoría, si se lleva al extremo, permite aproximar el comportamiento del aire al modelo fluidodinámico de Newton, usado con frecuencia para cálculos aproximados. -Capa entrópica Cuando se trata de cuerpos con nariz redondeada, además del efecto anterior, se da el de la existencia de una capa pegada al cuerpo de gran entropía. Esto se debe a la mezcla de dos cosas. Por una parte, el crecimiento de la entropía a través de la onda de choque, el cual es aún mayor cuanto más fuerte es la onda. Por otra parte, está el hecho de que cerca de la nariz, la onda de choque es prácticamente normal a la dirección de la corriente. Esto crea que un cierto número de líneas fluidas experimenten un gran crecimiento de la entropía mientras que otras no lo hacen, creando así una capa. Adicionalmente, esta gran entropía viene de la mano de una gran vorticidad según el teorema de Crocco. Esto supone una dificultad a la hora de definir las condiciones de contorno para el cálculo la capa límite. -Interacción viscosa Esto hace referencia a la interacción que existe entre el flujo no viscoso o ideal y el viscoso. El origen de esta interacción está en la energía cinética que posee el fluido a velocidades hipersónicas. Esta energía se transforma en energía interna al frenarse Introducción 3 el fluido en la capa límite, lo que se traduce en un aumento de la temperatura que, a través de la ley de Shuterland, se convierte en un aumento de la viscosidad del fluido como se puede ver en la Figura 1-1. Por último, este aumento en la viscosidad supone un aumento del grosor de la capa límite, es decir, de la zona viscosa. Cuanto más alto sea 𝑀∞, mayor será el grosor de la capa límite. Esto provoca que la forma del sólido virtualmente cambie provocando un cambio en la corriente no viscosa que a su vez cambia la capa límite, que a su vez, etc. Esta continua realimentación es lo que se conoce como la interacción viscosa. Figura 1-1 Ley de Shuterland. -Flujos de alta temperatura En relación con el apartado anterior, la extrema disipación de energía supone como ya se ha comentado un gran aumento de la temperatura cerca de la pared y esto puede desencadenar la disociación del aire o incluso su ionización. Además, dependiendo del material con el que esté fabricada la pared, esta también puede reaccionar. En resumen, se tiene una capa límite reaccionando químicamente. Esto afectará a la resistencia y la sustentación además de a las comunicaciones por radiación debido a que los electrones sueltos absorben la radiación de las ondas de radio. Además, para poder tener en cuenta estos efectos con precisión será necesario dejar a un lado la hipótesis de gas ideal. Como se ve en la Figura 1-2, a mayor temperatura, la diferencia entre ambos modelos se hace más acusada. Un cálculo preciso de la temperatura es en realidad crucial para asegurar el buen funcionamiento de todos los sistemas además de la supervivencia de las personas 4 Estudio de la Capas Límite en Régimen Hipersónico que haya dentro del vehículo, si las hubiese. Figura 1-2 Comparación de temperatura tras onda normal de choque según velocidad de vuelo y modelo de gas. Por último, las altas temperaturas provocan que, además de por convección, se transmita calor por radiación. Para que se tenga una referencia de la importancia que puede llegar a tener, en la reentrada de la misión Apollo supuso un 30% del calor total transferido. -Flujo de baja densidad La mayoría de los vuelos supersónicos se producen a muy grandes alturas, donde la atmósfera se torna muy débil. Este hecho provoca que la suposición del aire como un medio continuo deje de ser válida pudiendo provocar desde que tengan que ser revisadas las condiciones de contorno en la pared a la hora de calcular la capa límite hasta que las ecuaciones dejen del ser válidas por completo y sea necesario utilizar la teoría cinética para poder resolver el problema. Si bien este efecto es totalmente independiente del número de Mach, se ha creído interesante incluirlo en esta lista puesto que cualquier vehículo espacial en su Introducción 5 reentrada va a experimentarlo. Todos estos efectos pueden verse recapitulados en la Figura 1-3 proporcionada por [1]. De todos estos efectos, dado que este es el primer trabajo realizado en esta Escuela sobre este tema, solo algunos se van a considerar. Esencialmente va a realizarse un estudio de los efectos de compresibilidad en la capa límite, dejando para estudios posteriores los efectos de disociación, reacciones químicas y radiación. Figura 1-3 Efectos característicos del régimen hipersónico. 1.2. Relevancia Como se ha dicho en el apartado anterior, en el caso del régimen hipersónico, no existe una transición marcada como en paso de subsónico a supersónico sino que paulatinamente van apareciendo ciertos fenómenos físicos. Por tanto, cabe preguntarse si realmente es necesario tratar el régimen hipersónico como algo totalmente diferente. Para responder de forma rápida a esto pueden compararse los resultados obtenidos para un mismo estudio en régimen hipersónico y en supersónico. Lo más fácil es 6 Estudio de la Capas Límite en Régimen Hipersónico mirar la geometría de aviones (al final un avión es un objeto diseñado para sustentar) diseñados para volar en cada uno de estos regímenes y ver si existen diferencias sustanciales como existen en el caso de aviones subsónicos y supersónicos. Si observamos un avión supersónico cualquiera, sea militar como el Eurofigther Thyphoon o el Lockheed Martin F-35 Lightning II o civil como el Aérospatiale- BAC Concorde, claramente se sigue una estructura básica: un fuselaje aproximadamente cilidríco de sección circular con morro puntiagudo y alas en forma de ala delta de pequeño alargamiento. Figura 1-4 A la izquierda, fotografía del Concorde por Adrian Pingstone. A la derecha Eurofither Typhoon. Copyright Jeoffrey Lee. En cambio, para observar aviones hipersónicos hay que conformarse con proyectos aun experimentales como el X43A de la NASA o diseños que se quedaron en el túnel de viento como el Boeing X20-Dynasoar, también de la NASA. Figura 1-5 A la izquierda imagen conceptual del X-20 en los años 60 y a la derecha X43A en vuelo a principios de los 2000. Introducción 7 A pesar de haber más de 40 años entre ambos diseños donde se hicieron grandes avances en el campo del régimen hipersónico, hay detalles que siguen viéndose en ambos. El morro, al contrario que antes, ya no es puntiagudo sino redondeado llegando a poderse considerar romo. Además, todas las partes del avión (alas, fuselaje y motor) se integran en uno, manteniendo una superficie totalmente plana (la inferior en el X20 y la superior en el X43A). Esto último se observa mucho mejor en el segundo caso donde la aeronave es un cuerpo con dos dispositivos estabilizadores. Además de estos, otros muchos vehículos han sido diseñados para llegar y soportar al régimen hipersónico con el objetivo de ir más allá de la atmósfera y llegar al espacio. En esta tesitura, se pueden destacar dos tipos fácilmente reconocibles: el space shuttle y las cápsulas Apollo. Figura 1-6 Challenger despegando. Fotografía de la NASA. En el Challenger, que puede recordar al Dyna-Soar, pueden de nuevo verse las características destacadas antes: morro romo y panza totalmente lisa (una vez guarda el tren de aterrizaje). Por tanto, se hace patente la gran diferencia que existe entre en los dos regímenes puesto que, si no hubiese apenas cambios, los diseños tampoco serían tan diferentes. Es por ello que el régimen hipersónico debe considerarse y estudiarse como algo aparte del régimen supersónico. Pero esto no queda aquí. Todavía hay más razones que hacen necesario el estudio. La primera de ellas y en mi opinión la más importante actualmente es la reentrada en objetos desde el espacio a la Tierra ya sea con personas o no en su interior. 14 Estudio de la Capas Límite en Régimen Hipersónico Figura 1-8 Evolución del punto de desprendimiento de una capa límite. La separación de la corriente depende fundamentalmente de las condiciones exteriores siendo promovida por corrientes en los que existen gradientes adversos de presiones. La corriente tras el desprendimiento, sin apenas energía cinética, causa un gran cambio en el campo de presiones alrededor del sólido. Esto genera una resistencia conocida como resistencia de presión o de forma (debido a que el campo de presiones está íntimamente relacionado con la geometría del cuerpo). Esta resistencia es tanto mayor cuanto mayor sea la superficie donde la corriente está desprendida llegando al punto de que si la superficie es lo suficientemente grande, puede llegar a ser varios órdenes de magnitud mayor que la causada por la fricción. Figura 1-9 Distribución de presiones en cilindro circular para distintas corrientes. Introducción 15 En la Figura 1-9 se hace patente la importancia y la diferencia que supone solamente en el cálculo de la presión exterior el tener en cuenta la capa límite. Por tanto, la capa límite nos lleva a una doble fuente de resistencia: la de fricción en cualquier caso y el aumento de la de forma en el caso de que la corriente se vea desprendida. El desprendimiento por tanto es un efecto nada deseable y que en la aeronáutica intenta retrasarse lo más posible o que ni si quiera aparezca en ninguno de los componentes de una aeronave. Este punto de desprendimiento, como se verá más tarde, se adelanta cuanto más romo sea un cuerpo. Por otra parte, se dice que un cuerpo es romo cuando su longitud característica transversal a la corriente incidente es del mismo orden que la longitud paralela a la corriente. Por el contrario, se dice que un cuerpo es fuselado cuando la longitud transversal es mucho menor. Es por esto por lo que los perfiles aerodinámicos son tan esbeltos, aunque, como se ve en la imagen, a grandes ángulos de ataque, se comporta como un cuerpo romo. 1.4. Estructura del trabajo En el trabajo, antes de estudiar la capa límite, se estudiará la corriente exterior, necesaria para conocer las condiciones de contorno del tema central del trabajo. El cálculo de la corriente exterior de forma exacta es muy complicado. Gran parte de su complicación reside en el cambio de régimen que puede producirse tras la onda de choque, de subsónico a supersónico. Este cambio lleva consigo una alteración de las ecuaciones que gobiernan las ecuaciones, pasando de ser elípticas a ser hiperbólicas. Uno de los métodos que puede usarse es el método de las características, que lleva consigo el inconveniente de que sólo es útil en la zona supersónica del flujo. Además, para uso se requiere conocer los datos de los primeros puntos, lo que lleva a aproximarlos. Si el cuerpo de estudio es puntiagudo, toda la región fluida tras a onda será supersónico. Sin embargo, si la nariz es redondeada, antes de la zona supersónica existirá una subsónica, la cual non puede resolverse mediante el método de las características. 16 Estudio de la Capas Límite en Régimen Hipersónico Con todo esto, el método no es exacto debido a que es necesaria una resolución numérica que conlleva unos errores intrísecos. Por estas razones, en el presente trabajo se hará un cálculo simplificado de la corriente exterior. Se usarán métodos basados en la geometría del cuerpo a estudiar. Además, se presentará una modificación de los métodos presentes en la literatura básado en estos. Con esta modificación se pretende obtener mejores resultados sin que la complejidad aumente. El estudio de la capa límite también se hará de forma simplificada mediante suposiciones como que es muy fina y que no existe un gradiente de presiones en las líneas normales a la superficie. Se propone un método sencillo basado en el método de diferencias finitas donde se tendrá en cuenta la compresibilidad del aire. Por otra parte, no se tendrán en cuenta los efectos de disociación, reacciones químicas o radiación que pudieran ocurrir en el aire debido a las temperaturas que se alcanzan. Estos posibles efectos están recogidos en la Figura 1-10. Esto se dejará para futuros trabajos. Figura 1-10 Reacciones en el aire a 1 atmósfera. Equation Section (Next) Corriente exterior 17 2 CORRIENTE EXTERIOR omo ya se ha visto, la corriente exterior, a pesar de ser teóricamente más simple de estudiar debido a simplificaciones como la consideración de fluido ideal (μ=0) puede llegar a ser un estudio muy tedioso. Gran parte de la dificultad reside en que tras las ondas de choque, el aire puede pasar de estar en régimen subsónico a supersónico, lo que implica que las ecuaciones pasan de ser parabólicas a hiperbólicas, introduciendo así una dificultad extra. Figura 2-1 Esquema general de cuerpo en régimen hipersónico. Para evitar esto, se ha optado por utilizar métodos aproximados en el cálculo de la misma. No obstante, en este capítulo se seguirá el orden lógico: primero se tratarán las ondas de choque que permitirá conocer las condiciones iniciales de la corriente exterior y posteriormente se comentarán los distintos métodos que pueden usarse para estimar la corriente exterior alrededor de la superficie terminando por exponer cual es el método elegido para posteriormente estudiar la capa límite. 2.1. Ondas de choque Este apartado es quizás el más fácil de los problemas envueltos en el estudio del flujo alrededor del cuerpo. Con este estudio se pretende, a partir de las condiciones del aire en el infinito y la geometría del cuerpo objetivo de estudio obtener C 18 Estudio de la Capas Límite en Régimen Hipersónico relaciones que nos lleven a las condiciones justo detrás de la onda de choque. Dichas relaciones que se buscan son las ya conocidas relaciones de Rankine- Hugoniot (1-7 a 1-10) obtenidas a partir de balances de masa, cantidad de movimiento y energía asociados a un volumen de control infinitesimal en cada punto de la onda de choque que rodea y contiene a la misma. Como es conocido, tras cierto desarrollo algebeaico (ver [3]) se obtienen de a partir de ellas los saltos de las magnitudes a través de la onda como función del Mach incidente, M1, y del ángulo de incidencia 𝛽. Figura 2-2 Esquema general de onda de choque oblicua En particular, para la presión y la densidad se obtienen 22 21 1 ( 1) 1 2 1 sin pM p      , (2-1) 22 21 22 11 sin( 1) ( 1) 2 M M sin         (2-2) y haciendo uso de la ecuación de estado g PRT   , (2-3) Corriente exterior 19 se obtiene para la temperatura que 2 2 1 1 2 1 () () T p p T   . (2-4) Una vez conocidas las variables termodinámicas, puede obtenerse el salto de velocidades. Si u determina la velocidad según el eje x absoluto y v según el eje y, se tienen 22 21 2 11 2( 1 1) 1) ( sinuM Mu      , (2-5) 22 21 2 11 2( )cos 1( 1 ) sin M vM u      , (2-6) el valor del número de Mach 22 222 1 12 sin( si 1() 1 ))n( 2 sin M M M            , (2-7) así como la relación entre la presión de remanso tras la onda 𝑃𝑂2 y 𝑃1 1 2 22 1 11 12 1 o PP M PP          . (2-8) Adicionalmente, de las ecuaciones anteriores puede obtenerse otra que relaciona 𝑀1 y 𝛽 con el ángulo de deflexión de la corriente, 𝜃, : 22 1 2 1 1 2cot ( 2 ) t2 an M sin M cos         . (2-9) En la representación de la Figura 2-3 se observa que para cada número de Mach existe un ángulo de deflexión de la corriente (θ) máximo para que una onda oblicua adherida sea suficiente para la adaptación de la corriente. En el límite de número de Mach infinito se obtiene que es aproximadamente de 45º. 20 Estudio de la Capas Límite en Régimen Hipersónico Asimismo, se observa que para un valor dado de 𝜃>0 existen dos valores posibles de 𝛽, uno correspondiente a 𝑀2<0 (choque fuerte) y otro correspondiente a 𝑀2>0 (choque débil), siendo este último el que se observa en la naturaleza. Cuando 𝜃= 0, sin embargo, esta regla no se cumple y la onda de choque será normal, aunque la resolución de la ecuación pueda llevar a otro valor. Figura 2-3 Diagrama θ−β−M de [1]. Obsérvese que si 𝑀1≫1, se tiene el comportamiento asintótico 𝛽→1.2𝜃. Lo que significa que en el régimen hipersónico la onda de choque se ajusta bien (con un error del 20%) la geometría del cuerpo en la zona frontal. Corriente exterior 21 Figura 2-4 Bala y esfera en régimen hipersónico de [4]. 2.2. Métodos aproximados En este trabajo ya se ha comentado de la complejidad de resolución de la corriente exterior en régimen hipersónico. La resolución de las ecuaciones de forma analítica es imposible (solo podría llevarse a cabo haciendo simplificaciones sobre las mismas), lo que lleva a los métodos numéricos, que conllevan un error siempre. 22 Estudio de la Capas Límite en Régimen Hipersónico Para resolver las ecuaciones exactas numéricamente es necesario recurrir a métodos como el de las características. Estos métodos son muy complejos, siendo necesario suponer datos para poder empezar a usarlo. Además, debido al cambio de régimen existente en el caso de cuerpos redondeados (tras la onda el aire evoluciona desde subsónico a supersónico), estos métodos no son útiles en toda la zona. Esto lleva a que en este trabajo se usen métodos aproximados. Por simplicidad, y dado que como se verá más tarde los resultados son más que aceptables, se usarán diversos métodos basados en la geometría del objeto usando distintos mecanismos de expansión, unos basados en leyes newtonianas y otros en ondas de expansión. Método de Newton Este método fue propuesto por Newton para calcular la fuerza sobre una placa inclinada a baja velocidad. Más tarde se comprobó que era poco precisa pero irónicamente, tiene aplicación sobre flujos hipersónicos. El modelo de Newton consistía en que el fluido, modelado como una línea de corriente recta, cambiaba su velocidad normal a la superficie manteniendo la tangencial. De esta forma, si se tiene una placa inclinada un ángulo θ: Figura 2-5 Esquema teoría de impacto de Newton. Cambio en velocidad normal: 𝑉∞𝑠𝑖𝑛𝜃. Flujo de masa que incide en una superficie A: 𝜌∞𝑉∞𝐴𝑠𝑖𝑛𝜃. Tasa de cambio de cantidad de movimiento: 𝜌∞𝑉∞2𝐴𝑠𝑖𝑛2𝜃. Según la segunda ley de Newton, la tasa de cambio es igual a la fuerza sufrida por la placa. Esta fuerza no incluye el movimiento aleatorio de las partículas que originan p∞. Por tanto, 𝑝−𝑝∞=𝜌∞𝑉∞2𝑠𝑖𝑛2𝜃. Adimensionalizando, se tiene la ley de seno cuadrado para el coeficiente de Corriente exterior 23 presiones es 2 22sin 0.5 P PP CV       . (2-10) Según este método, el coeficiente de presiones sobre una superficie cualquiera 2D no depende de las condiciones exteriores del problema sino únicamente de la inclinación de la superficie en cada punto. Todos los métodos basados en la ley de Newton expuesta tienen como ventaja que es posible usarlos para cualquier tipo de cuerpo, ya sea redondeado o afilado. Sin embargo, presentan varias desventajas como que su precisión aumenta cuanto mayor sea el Mach de vuelo, provocando que la onda de choque se adhiera más al cuerpo, por lo que para bajos números de Mach tiene mayores errores. Por otra parte, dada la naturaleza balística del método solo puede usarse si 𝜃>0. Esto implica que la zona de sombra (𝜃<0), donde no impacta el fluido directamente se tiene una presión igual al último punto del cuerpo que no estaba en la zona de sombra. Ley de Newton modificada. En 1955, Lester Lees [5] propuso una alternativa al método de Newton. Según la ley de Newton, el coeficinete de presiones en el punto de remanso (𝜃=90º) tiene un valor de 2. Lester pensó que el método sería más acertado si en lugar de 2, tuviera el valor real al que se denominará 𝐶𝑝𝑚𝑎𝑥, quedando por tanto la ley de la siguiente forma: 2 sin max PP CC   . (2-11) Figura 2-6 Esequema general situación onda-cuerpo. 30 Estudio de la Capas Límite en Régimen Hipersónico mucho más pobres en supersónico. Esto es así por la simplicidad del modelo, la geometría de la onda de choque y el comportamiento de las ondas de expansión. Por una parte, existe el hecho de que las ondas de expasión que se crean en la superficiese propagan hasta la onda de choque, donde se ven reflejadas y vuelven al cuerpo. Si pintamos este dibujo para el caso supersónico y el hipersónico se tiene lo siguiente: Figura 2-13 Comparación de geometría de ondas en régimen supersónico e hipersónico. En el método no se contempla en ningún momento que las ondas se vean reflejadas y, como se ve en b, en régimen hipersónico, al estar la onda tan pegada al cuerpo, las ondas reflejadas aparecen al final del cuerpo mientras que en régimen supersónico están distribuidas por toda la superficie. Método híbrido Como se ha visto, cada uno de los métodos expuestos tiene sus ventajas y sus desventajas como el tipo de cuerpos al que pueden ser aplicados o la precisión conseguida. En un intento de aunar las ventajas de los diferentes métodos y suplir las deficiencias se propone un método diferente el cual recoge partes del método mejorado de Newton y el método de las ondas de expansión. El método supone un cuerpo con nariz redondeada en el que el método de las ondas de expansión no es aplicable debido a que el Mach tras la onda está por debajo de 1. Sin embargo, debido a la expansión, la velocidad aumenta llegando al punto en el que el régimen pasa a ser supersónico y es posible utilizar el método de las ondas de expansión a partir de ese punto. Corriente exterior 31 Las ventajas de este método es que la presión en el punto de remanso es extacta a la real y que a partir del punto en el que el se pasa a régimen supersónico la evolución se acerca más a la real. Sin embargo, no todos son ventajas. El principal de los inconvenientes del método es que, si se punto donde se produce el cambio de régimen difiere algo entre la realidad y la aproximación, en el método de las ondas de expansión aparece un error acumulado. Dado que el punto de remanso el valor de la presión es exacto, cuanto más cerca esté el punto de cambio de régimen, mejor será la aproximación puesto que el error acumulado hasta este punto será menor. Es por ello que se hace un estudio sobre la posición del punto donde Mach=1 en perfiles elípticos frente a dos parámetros: la relación de aspecto entre los semiejes y el Mach de vuelo. Figura 2-14 Dependencia del punto de cambio de régimen con los diferentes parámetros con 𝜎=0.2 a la izquierda y 𝑀∞=15 a la derecha. En la Figura 2-14 aparece representada la tendencia del punto. Se observa para empezar que al aumentar el Mach (lo que equivale a ser más hipersónico) el punto se mantiene por debajo del 1% de la cuerda del perfil y constante, haciéndose una vez más patente la independencia de la corriente exterior con el número de Mach cuando este es suficientemente grande. Por otra parte, el punto se retrasa cuanto menos esbelto es, llegando a estar al 5% de la cuerda, que, aunque es peor con respecto a otras relaciones de aspecto, sigue siendo un buen resultado. Como conclusión de este pequeño estudio, se obtiene que a priori se tendrán mejores resultados cuanto más esbelto sea el cuerpo y mayor el número de Mach. A modo de ejemplo, en la Figura 2-15 se presenta el coeficiente de presiones calculado por los dos métodos. Se ve que la diferencia no es muy grande, pero esa pequeña diferencia puede suponer un salto en la precisión. Aun así, sin datos experimentales no es posible corroborar que la precisión sea mejor al juntar ambos 32 Estudio de la Capas Límite en Régimen Hipersónico métodos para cuerpos de nariz redondeada. Figura 2-15 Comparación de coeficientes de presiones para el método de Newton y el método híbrido. 2.3. Cálculo del resto de variables termodinámicas, velocidad y resistencia a partir del campo de presiones En el anterior punto se han expuesto diferentes métodos con los que calcular el campo de presiones alrededor de un cuerpo sobre el que posteriormente se hará un estudio de la capa límite. Sin embargo, para poder hacer este estudio no basta con conocer el campo de presiones, es necesario tener todas las variables termodinámicas y la velocidad. Para calcularlas basta con saber que la presión que se conoce es a través de una línea de corriente, donde se conservan las magnitudes de remanso así como la entropía (es decir, es un movimiento isentrópico). Por tanto, basta con las ecuaciones de relaciones isentrópicas para obtener todo lo deseado. Dichas relaciones proporcionan 1)( / 21 1 o e e p Mp              , (2-17) Corriente exterior 33 2 11 2 oe e TM T    , (2-18) e e ge p TR   , (2-19) e e g e TMU R   , (2-20) Por último, va a ser de gran interés calcular la resistencia originada por la distribución de presiones sobre el cuerpo. Para ello, se integra la proyección sobre el eje X del campo de presiones: 0)(2 f x pX D P P dxn   . (2-21) Figura 2-16 Esquema de cuerpo para cálculo de resistencia de presión. Aun así, es más interesante el coeficiente de resistencia, puesto que es de mayor utilidad para estudios futuros. Para ello, basta con adimensinalizar la ecuación, quedando 0 2f DXp p f c C C n d l    . (2-22) 34 Estudio de la Capas Límite en Régimen Hipersónico 2.4. Distribución de presiones y resistencia de presión en distintos cuerpos bidimensionales. Antes de pasar al estudio de la capa límite recogido en el capítulo 3 es interesante observar el comportamiento de la corriente exterior en función tanto de la esbeltez como del Mach de vuelo. De la misma forma, también cobra una gran importancia calcular la resistencia aerodinámica que la distribución de presiones genera a lo largo del cuerpo. Una vez hecho el estudio en cada cuerpo se compararán los resultados entre ambos cuerpos. Los cuerpos que se estudiarán en este punto serán los mismos que se estudiarán en el capítulo 3: perfiles elípticos y parabólicos. Además, de forma separada se hará lo mismo con un cilindro y con una placa plana con borde de ataque semicircular. Perfil elíptico Para caracterizar la esbeltez de los perfiles elípticos se hace uso del parámetro 𝜎, que es la relación que existe entre el semieje menor y el mayor (b/a). De esta forma, 𝜎=0 correspondería a una placa plana y 𝜎=1 a un cilindro circular. Por otra parte, el barrido que se hará en Mach de vuelo irá desde 5 hasta 35. Figura 2-17 Esquema de perfil elíptico. La ecuación que describe la elipse (con a=1 para que sea adimensional) es 2 2 21 Y X   . Corriente exterior 35 El campo de presiones se calculará usando el método híbrido presentado en 2.2.5 ya que como se comentó, una vez el aire pasa a ser supersónico, existe un mayor fundamento físico en el proceso de expansión por medio de ondas de Pandtl-Meyer locales. Si se fija la relación de aspecto en 0.2 y se hace un barrido en Mach se obtiene la Figura 2-18, la cual puede tomarse como representativa del campo de presiones sobre una elipse en régimen hipersónico. En ella se observa primeramente que el gradiente de presiones sobre la superficie es siempre positivo, siendo muy acusado cerca del borde de ataque y suavizándose para manternerse prácticamente constante. Esto tiene grandes implicaciones en la capa límite ya que un gradiente siempre positivo significa que la corriente no se desprenderá. Además, a pesar de haber representadas 16 gráficas diferentes comprendidas entre Mach 5 y 35, apenas se ve diferencia entre ellas. Esto pone de manifiesto una vez más la independencia con el número de Mach que parece tener la corriente exterior como ya se vió cuando se presentó el método de Newton modificado. Figura 2-18 Coeficiente de presiones sobre elipse para distintos números de Mach con 𝜎=0.2. 36 Estudio de la Capas Límite en Régimen Hipersónico Si ahora en lugar de fijar la relación de aspecto se fija el número de Mach en 20, se obtiene la Figura 2-19 que permite analizar la influencia que tiene la relación de aspecto. Figura 2-19 Coeficiente de presiones sobre elipse a Mach 20 y diferentenes relaciones de aspecto. En primer lugar, se observa que el gradiente de presiones en el borde de ataque se suaviza con la relación de aspecto, debido a que la expansión necesaria es menor. Así mismo, esta tendencia cambia cuando se avanza en la cuerda llegando a ser mayor la presión cuando la relación de aspecto es pequeña a partir de aproximadamente la mitad de la cuerda. Por último, la transición de subsónico a supersónico puede apreciarse mucho mejor cuando 𝜎 es grande ya que el punto se retrasa y queda representado como un pico por el cambio de de mecanismo de expansión existente. Corriente exterior 37 La resistencia que genera las distribuciones de presiones siguen la misma tendencia como puede verse en la Figura 2-20. Figura 2-20 Coeficiente de resistencia de presión frente al número de Mach para diferentes 𝜎. Primeramente, como era de esperar, la resistencia aumenta con la relación de aspecto llegando aumentar su orden de magnitud en 2 hasta llegar a la unidad. Por otra parte, se vuelve a observar la independencia con el número de Mach a partir de Mach 8. La independencia con el número de Mach, a falta de futuros estudios con más detalle, tiene grandes implicaciones. Si tras hacer cálculos más rigurosos esta tendencia se mantuviese, implicaría que a la hora de hacer experimentos en túnel de viento bastaría con hacer los ensayos para un único número de Mach para obtener resultados tanto de variables termodinámicas como de velocidad y resistencia de presión para un amplio rango de velocidades hipersónicas. 38 Estudio de la Capas Límite en Régimen Hipersónico Perfil parabólico Para un perfil parabólico puede hacerse un estudio del mismo tipo y alcance que para los perfiles elípticos. Estos perfiles vienen descritos mediante la expresión 4 (1 ) max Y y X X    . En esta ocasión, el paramétro 𝜎 corresponderá al grosor máximo del perfil dividido por la cuerda (2ymax/c) siendo así análogo al parámetro usado para los perfiles elípticos. Figura 2-21 Esquema de perfil parabólico. Para este caso, a diferencia del anterior, no se usará siempre el método híbrido para calcular la corriente exterior. En lugar de ello, para casos los suficientemente esbeltos se usará el método de las ondas de expansión debido a los magníficos resultados que devuelve como se vió en la Figura 2-12. Cuando la relación de aspecto no permita usar este método, limitado por el ángulo del borde de ataque, pasará de nuevo a usarse el método híbrido. Los resultados de nuevo se presentarán de la misma forma que en el subpunto anterior para su estudio y posterior comparación. Corriente exterior 39 Figura 2-22 Coeficiente de presiones para perfil parabólico con 𝜎=0.2. Empezando con la dependencia con el número de Mach, en las figuras 2-12 y 2-13 se ha representado el coeficiente de presiones para dos relaciones de aspecto que requieren el uso de métodos diferentes. A pesar de la diferencia de métodos usados ambas mantienen dos propiedades en común. En primer lugar, a pesar de la diferencia en la evolución, ambas presentan un gradiente de presiones positivo a lo largo de todo el perfil (las diferencias entre ambos gradientes se verán más adelante) al igual que el presentado en los perfiles elípticos. Además, vuelve a poder observarse una tendencia a la independencia con el número de Mach cuando este es elevado. Esto puede verse observando como las líneas que corresponden a cada Mach se acercan cada vez más hasta prácticamente ser una. Esta tendencia se hace más acusada en el caso en el que se usa el método híbrido. 46 Estudio de la Capas Límite en Régimen Hipersónico Por otra parte, para el cono será necesario tener en cuenta la resistencia de base, es de decir la contribución que hace la base del cono. Teniendo en cuenta que es una base circular y que el ángulo de expansión es mayor que el que permiten las ondas de Prandtl-Meyer, el valor de la presión en esta zona es 0. De [8] se tiene una representación gráfica con resultados experimentales del coeficiente de resistencia tanto para una esfera como para un cono (Figura 2-32) con la que se comprobará la validez de la aproximación hecha. Figura 2-32 Coeficiente de resistencia experimental de un cono y una esfera de [8]. El cilindro con morro cónico y el el cono ofrecerán la misma resistencia debido a que la pared del cilindro no ofrece resistencia alguna dada su orientación con respecto a la corriente incidente. En la Figura 2-33 se observa que, a bajos números de Mach los resultados no son buenos, pero a partir de Mach 6 en la esfera, y Mach 8 en el cono, los resultados empiezan a concordar bastante bien. Corriente exterior 47 Figura 2-33 Coeficiente de resistencia aproximado para esfera y cono Es necesario tener en cuenta que, tanto el método híbrido usado para la esfera, como el método de ondas de expansión usado con el cono, tienen resultados muy pobres a bajos números de Mach como ya se había comentado cuando se presentaron, pero estos mejoran con el aumento del Mach. Equation Section 3 48 Estudio de la Capas Límite en Régimen Hipersónico Capa Límite 49 3 CAPA LÍMITE l estudio de la capa límite, último paso dentro del recorrido del trabajo, es de gran importancia debido a la crucialidad de sus resultados en aplicaciones reales del problema como puede ser el diseño de un vehículo espacial. Por una parte, del estudio de la capa límite se obtiene la resistencia originada por la fricción entre el aire y el objeto de estudio. Esta es conveniente compararla con la ya obtenida resistencia de presión para así determinar la importancia relativa de cada una de ellas. Por otra parte, se encuentran los efectos de temperatura. Estos, por la experiencia que se tiene en vuelo hipersónico, son los más importantes a tener en cuenta. Debido a los altos números de Mach, la temperatura cerca de los cuerpos aumenta mucho, y por consiguiente la temperatura de la superficie del cuerpo. En este trabajo no van a tenerse en cuenta varios de estos efectos como pueden ser las reacciones químicas que pueden producirse o la radiación. Sin embargo, de forma sencilla, se buscan, en lugar de resultados cuantitativos, resultados cualitativos que permitan hacernos una idea de la dependencia de las diferentes magnitudes con parámetros como la esbeltez o el Mach de vuelo. En el presente capítulo se presentarán las ecuaciones utilizadas para calcular la capa límite, así como el método numérico que se usará para implementar dichas ecuaciones. Para terminar, se mostrarán los resultados para diferentes estudios hechos de la capa límite para diferentes geometrías y velocidades de vuelo. Antes de poder presentar las ecuaciones es necesario hacer una pequeña introducción a ellas. En la capa límite, como cualquier movimiento de un fluido, según el número de Reynolds (el cual mide la importancia de las fuerzas de inercia frente a las de viscosidad) las partículas pueden moverse en régimen laminar o turbulento. Un movimiento es laminar si las partículas se mueven de forma ordenada siguiendo líneas de corriente que forman capas por encima del objeto. Por otra parte, el flujo será turbulento si deja de ser ordenado para ser caótico siendo el movimiento de las partículas aleatorio (No se estudia la velocidad total sino la velocidad media). E 50 Estudio de la Capas Límite en Régimen Hipersónico Figura 3-1 Perfil de velocidad en capa límite. Para decidir si el flujo es laminar o turbulento existen varios criterios. Para una placa plana resultado conveniente por su simplicidad el método de Pohlhausen. Este dice que la capa límite será turbulenta si el número de Reynolds local asociado al grosor de la capa es mayor que 645. Debido a que el grosor de la capa no es algo fijo, sino que va aumentando conforme recorre el objeto, la capa límite puede ser en una zona laminar y en otra turbulenta existiendo una zona de transición que en este trabajo se considerará infinitesimal. En el presente trabajo solo se harán cálculos para régimen laminar. Aún así, el método y los programas están preparados para implementar la transición a la turbulencia. 3.1. Capa límite laminar Aunque como se ha visto en el capítulo 1, las ecuaciones de la capa límite son en teoría las ecuaciones de Navier-Stokes (1-1 a 1-5) al completo ya que en esta ocasión la viscosidad es importante por pequeña que sea, otro tipo de simplificaciones pueden hacerse. Estas están relacionadas con el grosor de la capa. Por comodidad del lector se reproducen aquí de nuevo las ecuaciones 1-1 a 1-5 con ciertas simplificaciones ya hechas: Capa Límite 51 ()0V t     , (3-1) 22 22 )( () p tx u u u uV xy                , (3-2) 22 22 (() )p v v ty vvV xy              , (3-3)   v eV e p V K T t               , (3-4) Donde se han despreciado los esfueros de fuerzas másicas al ser 𝑈∞ 2≫𝑔𝑙, donde l es la longitud caracterísitica del cuerpo. Pueden sustituirse los términos de fuerza viscosa pueden sustituirse por la ley de Newton, que los relaciona con el gradiante de velocidad como se comentó en la sección 1.3.3. Adicionalmente se ha eliminado la ecuación de movimiento proyectada en el eje z ya que como ya se ha comentado, en este trabajo solo va a tenerse en cuenta el problema bidimensional. Adicionalmente, es necesario usar la ecuación de estado, así como calcular los valores de 𝜇 y 𝑘. Se utilizará la ley de Shuterland para la viscosidad, quedando así las respectivas expresiones: wT S T T T S               , (3-5) eg PTR   , (3-6) Pr p kc   . (3-7) Una vez terminada la adecuación de las ecuaciones, se pasa a las simplificaciones de capa límite propiamente dichas. Debido a que los términos convectivos y viscosos deben ser del mismo orden, pueden compararse para averiguar el orden de magnitud del grosor de la capa límite en comparación con la longitud característica del problema. Para ello, las magnitudes físicas caracterñisticas son las siguientes: 𝑥~𝐿 𝑦~𝛿 𝑢~𝑈∞ 𝑣~𝑉𝑦𝑐 𝜌~𝜌∞ 52 Estudio de la Capas Límite en Régimen Hipersónico Comparando el termino convectivo en x y el viscoso en y se tiene que 2 1 ~ ~ 1 UU UL L LU Re        . (3-8) A partir de esta simplificación, otras pueden hacerse: teniendo en cuenta los órdenes de magnitud, es posible despreciar los términos viscosos en x en comparación con los mismos en y. Además, debido al pequeño grosor de la capa, las variaciones de la presión en la dirección y son mucho menores que en la dirección x. Por tanto, la ecuación 1-3 puede dejarse al un lado, siendo solo utilizada en ocasiones para hacer correcciones. Al no haber variaciones de la presión en el eje y, en cada estación la presión será igual a la de la corriente exterior. La ecuación 3-4 hace referencia a la capa límite térmica, que da información sobre el perfil de temperaturas dentro de la capa límite. Esta ecuación también puede simplificarse y expresarse de forma alternativa. Por una parte, la variación de energía se puede expresar como variación de entalpia, equivalente a cp multiplicada por la temperatura. Por otra parte, los términos convectivos deben ser del mismo orden que los de conducción. Para ello, es necesario que se de lo siguiente: 𝛿𝑇 𝑙~1√𝑅𝑒𝑃𝑟. Dado que Pr es de orden unidad para el aire, de nuevo el grosor de esta capa es mucho menor que la longitud característica del problema. Esto tiene como consecuencia que se pueden llevar a cabo simplificaciones de la misma índole que la ecuación de cantidad de movimiento. Por tanto, la ecuación de la energía queda de la siguiente forma: 2 e pp u dp T T T cc y u v u k x y dx y y                      . (3-9) Para la resolución de las ecuaciones se hacen necesarias condiciones de contorno en ambos extremos tanto en velocidad como en temperatura. En la pared se tiene que 𝑢=0 , 𝑣=0, 𝑇=𝑇𝑝 (𝑜 𝜕𝑇 𝜕𝑦=0), Capa Límite 53 según la pared sea o no adiabática. Por su parte, para el acoplamiento con la corriente exterior es necesario que se cumpla en el último punto que 𝑢=𝑢𝑒, 𝑇=𝑇𝑒. Por último, se adimensionalizan las ecuaciones. De esta forma se consigue que los resultados no dependan del tamaño del cuerpo, de la velocidad, altura o número de Reynolds al que vuele, solo de la forma, Mach y modelo de aire usado. Para ello, se introduce el escalado de varibles 𝑢→𝑈∞𝑢, 𝑣→ 𝑈∞𝑣, 𝑥→𝑙𝜉, 𝑦→ 𝑙𝜂 √𝑅𝑒 , 𝑢𝑒→𝑈∞𝑢𝑒, 𝑝→𝑝∞𝑝, 𝑇→𝑇∞𝜃, 𝜌→𝜌∞𝜌, 𝜇→𝜇∞𝜇, 𝑘→𝑘∞𝑘, donde las variables del segundo término son adimensionales Por consiguiente, las ecuaciones de la capa límite laminar adimensionalazas son )( )( 0 uv      , (3-10) 2 1e dp u u u uv Md                          , (3-11) 2 2 11) P( 1 r e dp u u v u k M d                                              (3-12) y adicionamelte, 𝑝𝑒 𝜌=𝜃, 𝑘=𝜇, 𝜇=𝜃𝑤(1+𝛿1 𝑇∞ 𝜃+𝛿1 𝑇∞)𝜎. 54 Estudio de la Capas Límite en Régimen Hipersónico 3.2. Resolución numérica de las ecuaciones: método de las diferencias finitas El método de las diferencias finitas es un método de resolución de ecuaciones diferenciales basado en la discretización de las derivadas en el recinto donde quieren resolverse las ecuaciones. Este método es general y adecuado para el sistema de ecuaciones bidimensionales que constituye el problema de la capa límite. Mallado del plano En primer lugar, es necesario discretizar el recinto en el que va a resolverse a través de una malla en los ejes x, que sigue la superficie del cuerpo y el eje y, perpendicular a la superficie. Con esta discretización se obtiene un sistema matricial de cuya resolución se realiza estación a estación, siendo una estación la línea de puntos vertical con un valor de x (denotado con el subíndice n) y todos los valores de y pertenecientes a esa línea (denotados con subíndice j). El resultado del problema será el campo de velocidades (u) y de temperaturas (𝜃) en todos los puntos de la malla. Debido a la resolución numérica de las ecuaciones, no es posible tener el campo de velocidades y temperaturas al completo. Para obtener el mayor número posibles de valores bastará con hacer una malla con más puntos, teniendo esto un coste de tiempo en la resolución. Por tanto, se tendrá que llegar a una solución de compromiso en la que la obtención de buenos resultados no conlleve un excesivo coste temporal. El mallado puede ser de dos formas: uniforme o variable. En el caso del mallado uniforme, los puntos de la malla son equiespaciados en ambos ejes (aunque el espaciado de cada eje no tiene que ser el mismo). Este es un mallado que es adecuado si no se desea obtener mejores resultados en una zona en especial sino tener la misma precisión en todo el espacio. Capa Límite 55 Figura 3-2 Mallado uniforme con L puntos en el eje x y M puntos en la dirección y. Por otra parte, un mallado variable es aquel en el que el alguno de los ejes (o ambos) no presenta líneas equiespaciadas. Esto es de gran utilidad si se desea o requiere una mayor precisión en alguna zona en concreto. Por ejemplo, en el caso de perfiles elípticos, es importante tener un mallado muy fino en los primeros puntos ya que el borde redondeado causa que las variables cambien mucho en esa zona. Este mallado variable puede conseguirse haciendo que los puntos del eje no sigan una ley lineal sino cuadrática, cúbica… Cuanto mayor sea el grado, mayor será la diferencia de espaciado entre los primeros puntos y los últimos. Una posible desventaja de este tipo de mallado es, que si el número de puntos no es el suficiente (o el grado del polinomio es muy grande), puede suceder que la precisión en los últimos puntos se vea muy limitada. Figura 3-3 Mallado variable en x. 62 Estudio de la Capas Límite en Régimen Hipersónico Figura 3-7 Perfil de temperaturas en placa plana adiabática. Por último, para confirmar la veracidad de los datos puede compararse el coeficiente de resistencia de fricción obtenido con el reportado por Van Driest en [9]. 3.4. Resistencia de fricción La resistencia de fricción es la consecuencia de la condición de no deslizamiento de la pared. Según la ley de Newton, el tensor de esfuerzos en la pared es proporcional a la viscosidad y el gradiente de velocidad en la pared como se indicaba en la ecuación (1-15). A partir del mismo, el coeficiente de fricción en cada punto se tiene que es Capa Límite 63 0 2 Re f Cu       . (3-15) Para obtener el coeficiente de resistencia de fricción basta con integrar el coeficiente fricción múltiplicado por el vector tangente a la superficie proyectado en el eje X a lo largo de la superficie, teniéndose que 0 f D f f X C C t d   (3-16) Al igual que se ha hecho con los perfiles de velocidades y temperatura, es posible comparar los resultados de resistencia en la placa plana. En la Figura 3-8se observa los resultados presentados en [9]. Figura 3-8 Coeficiente de resistencia de fricción sobre placa plana laminar de [9] Por su parte, en la Figura 3-9 se tienen los resultados obtenidos por MATLAB. Existe una buena concordancia en la mayoría de los casos. Se observa que la tendencia de la resistencia es a bajar cuando el número de Mach crece, así como cuanto mayor sea la temperatura de la placa en relación con la exterior. 64 Estudio de la Capas Límite en Régimen Hipersónico Figura 3-9 Coeficiente de resistencia de fricción en una placa plana para diferentes relaciones de temperaturas. Sin embargo, los perfiles objeto de estudio no presentan la misma evolución que la placa plana. En el caso de una elipse se obtiene la Figura 3-10. En ella se observa que, aunque existe una variación con el número de Mach, esta es muy pequeña (se ve aumentada por la escala de la representación). Así mismo se observa que la resistencia disminuye cuanto más roma es la elipse ( 𝜎 crece). Esto se debe a que el gradiente de presiones originado por la expansión aumenta con la curvatura del cuerpo, provocando que, a bajas relaciones de aspecto, donde el cambio de curvatura es más acusado, el grosor de la capa disminuya y aumente el gradiente de velocidad en la pared. Capa Límite 65 Figura 3-10 Coeficiente de resistencia de fricción en perfiles elípticos A bajos números de Mach se observa un crecimiento de la resistencia con 𝑀∞. Esto se debe a que el gradiente favorable de presiones aumenta con el Mach, ocasionando que el gradiente de velocidad sea mayor en la pared. Esto puede comprobarse en la Figura 3-11, donde se han representado la evolución del perfil de velocidades a dos números de Mach distintos. Figura 3-11 Evolución de perfil de velocidades a Mach 5 (izqda) y 25 (drcha). 66 Estudio de la Capas Límite en Régimen Hipersónico La disminución de la resistencia a altos números de Mach está causada por el aumento del grosor de la capa, que convella una dismincuión del gradiende velocidad. Este aumento del grosor está asociado al aumento del calor generado por fricción. Esto puede verse en la Figura 3-12 donde se ha representado la evolución de la capa límite térmica. Figura 3-12 Evolución de perfil de temperaturas a Mach 5 (izqda) y 25 (drcha). En ella se ve que, aunque a simple viste ambas parecen iguales, existe un factor de escala muy importante debido al cambio de grosor y a la temperatura que se alcanza en la capa. Esto provoca que a altos números de Mach se tenga un mayor geradiente de temperaturas en la pared. Por su parte, los valores de la resistencia de fricción perfiles parabólicos vienen recogidos en la Figura 3-13. Cuantitativamente, se observa que en los perfiles parabólicos la resistencia es menor también en este caso (en la sección 2 se vió un comportamiento similiar para la resistencia de presión). Adicionalmente, la dependencia con el Mach es muy similiar, presentando el segundo caso su máximo a menores números de Mach. Capa Límite 67 Figura 3-13 Coeficiente de resistencia de fricción en perfiles parabólicos. La justificación de estos resultados no difiere de la dada para los perfiles elípticos. Figura 3-14 Evolución de perfil de velocidades a Mach 5 (izqda) y 25 (drcha) En la Figura 3-14 y la Figura 3-15 se observa la evolución tanto del perfil de velocidades como del de temperaturas para dos números de Mach en perfiles 68 Estudio de la Capas Límite en Régimen Hipersónico parabólicos con los que poder justificar los resultados obtenidos. Figura 3-15 Evolución de perfil de temperaturas a Mach 5 (izqda) y 25 (drcha). Por último, es interesante conocer el valor de esta resistencia en el caso de que la pared sea adiabática, es decir, no exista transferencia de calor. Este caso sería equivalente a mantener el vuelo el tiempo suficiente para que se llegase al equilibrio térmico en la superficie. En este caso, el coeficiente de fricciíon obtenido para perfiles elípticos es el que se observa en la Figura 3-16. Figura 3-16 Coeficiente de resistencia de fricción en perfiles elípticos de pared adiabática. Capa Límite 69 En ella se obsera que la dependencia con la relación de aspecto es creciente, aunque estre crecimiento disminuye mucho llegando a ser prácticamente constante. Por su parte, el aumento de Mach conlleva el crecimiento de la resistencia también. Por su parte, el caso de los perfiles parabílicos viene representado en la Figura 3-17. En ella se ve que la dependencia con el número de Mach sigue siendo creciente. Sin embargo, existen dos tramos en la dependencia con 𝜎. Este cambio de tendencia de creciente a decreciente se debe al cambio físico que sufre el problema: por encima de 𝜎=0.4 la onda deja de ser oblicua separanadose en la parte frontal (razón por la que el método de ondas de expansión deja de ser útil) y la evolución del aire se ve modificada. Figura 3-17 Coeficiente de resistencia de fricción para perfiles parabólicos de pared adiabática. A la vista de los resultados obtenidos, y comparando con los valores de la resistencia de presión, se observa que, para cuerpos muy esbeltos, ambos coeficientes tienen valores pequeños y del mismo orden de magnitud. Sin embargo, para cuerpos romos, la resistencia de fricción apenas sufre cambios mientras que la de presión se vuelve mucho más importante por llegar a tener orden unidad. Por tanto, las resistencias serán equiparables solo en el caso de que el cuerpo sea esbelto. 70 Estudio de la Capas Límite en Régimen Hipersónico 3.5. Transferencia de calor La transferencia de calor existente entre el aire y un cuerpo es una de las magnitudes críticas en el régimen hipersónico. Es bien sabido que en una de las aplicaciones actuales más comunes en la que se alcanza velocidades hipersónicas es en la reentrada de vehículos espaciales y que las altas temperaturas han sido causa de más de un accidente. Es por eso que se hace especial atención a este apartado dentro de la sección de capa límite. Existen varios mecanismos por los que el calor puede transferirse del aire al cuerpo o viceversa: conducción, convección y radiación. De los tres, en este trabajo, se va a despreciar el calor transferido por radiación por simplicidad del modelo. El calor total transferido por unidad de tiempo y envergadura (el transferido en un perfil con espesor tendiendo a 0) puede calcularse integrando a lo largo de la superficie el calor transferido en cada punto. De esta forma, en variables físicas, se tiene: 0 f xT Q K dx y        (3-17) En la ecuación, que puede adimensionalizarse de la misma forma que se hizo en la sección 3.1, x corresponde con la usada para el cálculo de la capa límite. Además de presentar el calor transferido por los cuerpos y la dependencia con el número de Mach y la relación de aspecto, también es interesante conocer la importancia del calor trasnferido mediante convección frente al transferido mediante conducción. Esto se consigue a través del número de Nusselt, que se define como el coeficiente de película (que aparece en el calor convectivo) partido por la conductividad térmica por unidad de longitud (K/c). / h Nu kc   (3-18) Para calcularlo, se supone que el calor calculado anteriormente es equivalente al de la ecuación (3-19), donde c corresponde a la cuerda del perfil. Igualando ambas, se puede obtener el número de Nusselt. El calor de convección es () p Q h c T T     . (3-19) El calor, en variables adimensionales viene expresado por Capa Límite 71 0 Re f Q k T k d        . (3-20) Si esto lo igualamos a (3-19), se tiene la igualdad 0 Re ( ) f p k T k d h c T T             . Por último, manipulando la igualdad se llega a la expresión del número de Nusselt deseada: 0 Re () f p T Nu k d TT         (3-21) La veracidad de los reultados que vuelca el programa puede comprobarse mediante la literatura. En ella pueden encontrarse dos casos: el de un cilindro y el de una placa plana con borde redondeado. De ahí que se estudiase también la corriente exterior sobre ellos. En el caso de un cilindro, los resultados reportados por la literatura y por el programa de MATLAB están plasmados en la Figura 3-18. En ellas se ve que se obtienen resultados muy acertados del programa. Lo único que difiere es que aparece una pequeña dependencia con el Mach y un pico correspondiente al cambio de método usado para calcular la corriente exterior. Figura 3-18 Distribución de transferencia de calor en cilndro (a la izquierda resultados numéricos, a la derecha resultados de la literatura). 78 Estudio de la Capas Límite en Régimen Hipersónico Al igual que con el coeficiente de resistencia de friccion, la dependencia con el Mach es creciente mientras que la dependencia con la relación de aspeto es creciente hasta el cambio de físico del problema, donde comineza a decrecer para llegar a ser prácticamente constante. Por otra parte, la trasnferencia de calor total está representada en la Figura 3-28. Puede observarse que, a parte de presentar la misma dependencia con los parámetros que en el caso anterior, su valor es considerablemente menor. Figura 3-28 Transferencia total de calor en perfiles parabólicos Por último, resta conocer la importancia de la transferencia de calor por ceonvección frente a la generada por la fricción. Esto viene representado en la Figura 3-29 a través del número de Nusselt. Capa Límite 79 Figura 3-29 Dependencia del número de Nusselt con los diferentes parámetros. La dependencia con los parámetros es muy similar a la obtenida para los perfiles elípticos. Sin embargo, su valor es menor. Esto indica que en perfiles parabólicos la transferencia por conducción es más importante, debido a que cerca del borde ataque la pendiente difiere mucho de los perfiles elípticos, y por consigiuente la corriente exterior y la capa límite. Resulta interesante, además de estudiar la transferencia de calor, estudiar a que temperatura tienden los diferentes puntos del perfil en el límite en el que se alcanza el equilibrio térmico entre el aire y el perfil. Esta situación es equivalente a considerar que la pared es adiabática. En perfiles elípticos se tienen que la temperatura en la pared decrece a lo largo del perfil teniendo un gran crecimiento con el número de Mach como puede verse en la Figura 3-30. 80 Estudio de la Capas Límite en Régimen Hipersónico Figura 3-30 Temperatura en perfil elíptico (𝜎=0.6) adiabático. En la Figura 3-31 se observa cómo evoluciona la temperatura para un mismo número de Mach y diferentes relaciones de aspecto. La razón por la que la temperatura en la parte frontal y la final no dependen de la relación de aspecto es que los métodos usados dependen de la pendiente de la superficie y en esas zonas la pendiente depende poco de la relación de aspecto. Figura 3-31 Temperatura en perfiles elípticos (M=21) adiabático. Capa Límite 81 Los resultados para perfiles parabólicos están recogidos en las figuras 3-32 y 3-33. En ellas puede verse que no existe gran diferencia con el caso de perfiles elípticos en cuanto a la evolución y los valores numéricos. Figura 3-32 Temperatura en la pared en perfile parabólico (𝜎=0.6) para diferentes números de Mach Figura 3-33 Temperatura en la pared de difernetes perfiles parabólicos (𝑀∞=21). Cabe recordar que de este estudio que se ha hecho sólo se debe tomar la evolución de los parámetros como válida debido a las simplificaciones hechas en la resolución del problema. 82 Estudio de la Capas Límite en Régimen Hipersónico Conclusiones 83 4 CONCLUSIONES n el presente trabajo se ha querido hacer un primer estudio del régimen hipersónico con vistas a servir de base y dar pie a su desarrollo en trabajos posteriores. Para estudiar la corriente exterior no viscosa se ha simplificado el problema usando métodos basados en la geometría del cuerpo en lugar del sistema de ecuaciones completo. Además, se ha omitido la posible interacción entre corriente exterior y la capa límite. De esta forma, el objetivo del estudio de la corriente exterior se reduce al de calcular la resistencia de presión, así como la evolución de las diferentes magnitudes y su dependencia con la esbeltez y el Mach. Además, los resultados han sido usados como condiciones de contorno para el estudio de la capa límite. Además, se ha propuesto un método de resolución alternativo hasta ahora desconocido en la literatura, bautizado como método híbrido, que permité usar el método de ondas de expansión fuera de sus limitaciones originales (borde de ataque redondeado y esbeltez mínima necesaria) con el objetivo de usar un método con mayor fundamento físico. Por último, se ha hecho un estudio de la capa límite mecánica y térmica, que están acopladas. Haciendo uso de métodos sencillos basados en el de diferencias finitas para resolver las ecuaciones de Navier-Stokes se han conseguido resultados con la precisión de métodos más complejos basados en la autosemejanza usados en la literatura. Efectos de alta temperatura como las reacciones químicas o la radiación han quedado fuera del estudio, estando los programas preparados para la implementación de las ecuaciones necesarias para que sean tenidos en cuenta. Dentro de la capa límite se han hecho estudios de la resistencia de fricción y la transferencia de calor. Los resultados de ambos deben ser interpretados de forma cualitativa y no cuantitativa debido a las simplificaciones hechas. E 84 Estudio de la Capas Límite en Régimen Hipersónico Anexos 85 ANEXOS A Códigos de Matlab A.1 Método de Newton mejorado function [Pe_Pinf,C_P]= Newton_mej_v2(Minf,gamma,N,theta) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%% C_Pmax %%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % C_Pmax = 2/(gamma*Minf^2)*(((gamma+1)^2*Minf^2/(4*gamma*Minf^2- 2*(gamma-1)))^(gamma/(gamma-1))*((1- gamma+2*gamma*Minf^2)/(gamma+1))-1); % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%% C_P %%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % for j=1:N if theta(j)>=0 C_P(j) = C_Pmax*sin(theta(j)).^2; m = j; else C_P(j) = C_P(m); end end % Pe_Pinf = C_P*0.5*gamma*Minf^2+1; % End 86 Estudio de la Capas Límite en Régimen Hipersónico A.2 Método de Ondas de choque function [Ue_Uinf,Pe_Pinf,Te_Tinf,Roe_Roinf]=ShockExp(Minf,Nx,model o_aire,theta) gam=1.4; M=linspace(1,50,1000); aux=sqrt(M.^2-1); aux1=sqrt((gam+1)/(gam-1)); f=aux1*atan(aux/aux1)-atan(aux); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% Condiciones exteriores %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if modelo_aire==0 %Gas ideal gamma = 1.4; Rg = 287; cp = 1004; elseif modelo_aire==1 %Alta temperatura %completar end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Onda de Choque (Cond. Ini.) %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% teta=theta(1); par=[Minf teta]; TOLa=0; TOLr=1e-11; maxit=50; x0=1.2*teta; [beta,iflag,iter,normf]=newtoncompleto(@ThetaBetaMach,x0,T OLr,TOLa,maxit,par); %%Resulve la ecuación que relaciona theta y beta if iflag==1 iflag=1 end [U2_U1,Mn2,P2_Pinf,T2_Tinf,Ro2_Roinf]=onda_choque(Minf,bet a,modelo_aire); Anexos 87 M2=Mn2/sin(beta-teta); %Necesario calculcar Mach y Presion tras la onda Mx(1)=M2; %%%%%%%%%%%%%%%%%%%%%%%%%% %%% Evolución del Mach %%% %%%%%%%%%%%%%%%%%%%%%%%%%% for i=2:Nx, M1=Mx(i-1); f1=interp1(M,f,M1); dtheta=theta(i-1)-theta(i); f2=f1+dtheta; M2=interp1(f,M,f2); Mx(i)=M2; end % plot(x,Mx) Pe_Pinf(1)=P2_Pinf; Pe_Pinf(2:Nx)=Pe_Pinf(1)*(1+0.2*Mx(1)^2)^(gamma/(gamma- 1))./(1+0.2*Mx(2:Nx).^2).^(gamma/(gamma-1)); m=1; while m<=Nx if Pe_Pinf(m)<0.00001 Pe_Pinf(m)=0; end m=m+1; end Te_Tinf=T2_Tinf*(1+(gamma-1)/2*Mx(1).^2)./(1+(gamma- 1)/2*Mx.^2); % Roe_Roinf= Pe_Pinf./Te_Tinf; % Ue_Uinf=Mx.*sqrt(Te_Tinf)/Minf; % end 94 Estudio de la Capas Límite en Régimen Hipersónico %TEMPERATURA % At=spdiags(urhonm1',0,Ny,Ny)/hx(n)+spdiags(vrhonm1'- 1/Pr*dcappaynm1',0,Ny,Ny)*Dy- 1/Pr*spdiags(cappanm1',0,Ny,Ny)*Dy2; %Tp impuesta % rt(1)=tp; % At(1,:)=0; At(1,1)=1; % Adiabatico: rt(1)=0; At(1,:)=0; At(1,1)=f1; At(1,2)=f2; At(1,3)=f3; % At(Ny,:)=0; At(Ny,Ny)=1; rt(2:Ny-1)=urhonm1(2:(Ny-1)).*tnm1(2:(Ny- 1))/hx(n)+Minf^2*(gamma-1)*munm1(2:(Ny-1)).*(Dy(2:Ny- 1,:)*un')'.^2+(gamma-1)*unm1(2:Ny-1).*dpe_dx(n)/gamma; %falta añadir el termino del gradiente de presiones rt(Ny)=te(n); tn(1:Ny)=(At\rt')'; figure(1) hold on plot(tn,yv,'b') pause(0.01) hold off % %TRANSFERENCIA DE CALOR dtdy_0(n)=(f1*tn(1)+f2*tn(2)+f3*tn(3))*sqrt(Re).*cappanm1( 1); q(n)=(f1*tn(1)+f2*tn(2)+f3*tn(3))*sqrt(Re); % %Densidad, viscosidad: rho(1:Ny)=pe(n)./tn(1:Ny); mumol(1:Ny)=tn.^omega.*((1+CSuth)./(tn+CSuth)).^sigma; % dudy_0=f1*un(1)+f2*un(2)+f3*un(3); Cf(n)=(2/sqrt(Re))*mumol(1)*dudy_0; % % Cálculo de vn durhon(1:Ny)=(un.*rho-urhonm1)/hx(n); vrhon(1)=0; for j=2:Ny, Anexos 95 vrhon(j)=vrhon(j-1)-0.5*(yv(j)-yv(j- 1))*(durhon(j)+durhon(j-1)); end vn=vrhon./rho; % a estación siguiente unm1=un; vnm1=vn; tnm1=tn; rhonm1=rho; mumolnm1=mumol; end %CALCULO DE CD CDf=trapz(xv,Cf.*tX); 96 Estudio de la Capas Límite en Régimen Hipersónico 11 REFERENCIAS [1] J. D. Anderson Jr., Hypersonic and High Temperature Gas Dynamics. Second Edition, AIAA, 2006. [2] J. D. Anderson Jr., «Ludwing Prandtl's Boundary Layer,» Physics Today, 2005. [3] A. Barrero Ripoll y M. Pérez-Saborid Sánchez-Pastor, Fundamentos y Aplicaciones de la Mecánica de Fluidos, Sevilla: McGraw Hill, 2005. [4] A. H. Shapiro, The dynamics and thermodynamics of compressible fluid flow, New York: The Ronald Press Company, 1954. [5] L. Lees, «Hypersonic Flow,» de Fifth International Aeronautical Conference, Los Angeles, New York, Inst. of Aeronautical Sciences, 1955, pp. 241-276. [6] R. N. Cox y L. F. Crabtree, Elements of Hypersonic Aerodynamics, New York: Academic Press, 1965. [7] A. J. Eggers , C. A. Syvertson y S. Kraus, «A Study of Inviscid Flow about Airfoils at High Supersonic Speeds,» NACA Report 1123, 1953. [8] Gran Enciclopedia Soviética, 1979. [9] E. R. Van Driest, «Invesstigation of Laminar Boundary Layer in compressible fluids using the Crocco Method,» National Advisory Committee for Aeronautics, Washington, 1952. [10] J. Peral González, Cálculo de la capa límite compresible sobre perfiles aerodinámicos, Sevilla: Escuela Técnica Superior de Ingenieros, 2015. [11] J. D. Anderson Jr., Fundamentals of Aerodynamics, New York: McGraw-Hill series in aeronautical and aerospace engineering, 2011. 12 12 Estudio de la Capas Límite en Régimen Hipersónico