Area de mecanica de fluidos universidad de oviedo

1. Universidad de OviedoÁREA DE MECÁNICA DE FLUIDOShttp://web.uniovi.es/Areas/Mecanica.Fluidos/ MONOGRAFÍAS DE MECÁNICA DE FLUIDOSTÉCNICAS NUMÉRICAS EN MECÁNICA DE FLUIDOS Rafael Ballesteros Tajadura José González PérezJesús Manuel Fernández OroKatia María Argüelles Díaz 2. © Rafael Ballesteros TajaduraISBN 84-607-9546-2DEPÓSITO LEGAL AS-05282-2003Gijón 2003 3. Técnicas Numéricas en Mecánica de Fluidos 3Índice.Prólogo .................................................................................... 51.- Introducción....................................................................... 7 1.1.- Perspectiva histórica .................................................................................................. 7 1.2.- Aplicaciones habituales de las técnicas numéricas ................................................... 8 1.3.- Niveles de aproximación empleados en las técnicas numéricas. ........................102.- Modelos físico-matemáticos............................................ 13 2.1.- Modelos de flujo potencial ..................................................................................13 2.2.- Ecuaciones para flujo ideal..................................................................................13 2.3.- Modelo para flujo ideal, estacionario y rotacional ..............................................14 2.4.- Solución directa de las ecuaciones de Navier-Stokes..........................................14 2.5.- Modelos parabólicos de las ecuaciones de Navier-Stokes ..................................14 2.6.- Modelo de flujo incompresible............................................................................15 2.7.- Modelos para la simulación de la turbulencia .....................................................15 2.8.- Modelos de dos ecuaciones: modelo k-ε .............................................................18 2.9.- Modelos matemáticos para la Capa Límite .........................................................203.- Aspectos matemáticos del procedimiento de cálculo ...... 23 3.1.- Condiciones iniciales y de contorno....................................................................25 3.2.- Consideraciones físicas .......................................................................................274.- Métodos de resolución y naturaleza de las ecuaciones degobierno ................................................................................ 29 4.1.- Ecuaciones de gobierno, condiciones de contorno y condiciones iniciales ........30 4.2.- Discretización de las ecuaciones .........................................................................31 4.3.- Métodos de discretización empleados .................................................................38 4.4.- Discretización temporal.......................................................................................56 4.5.- Tipos de aproximación numérica utilizadas en la práctica..................................58 4.6.- Resolución numérica de problemas sencillos......................................................62 4. 4Técnicas Numéricas en Mecánica de Fluidos5.- Discretización del dominio: generación de mallados ...... 99 5.1.- Clasificación de los mallados basada en la conectividad y estructura de datos......................................................................................... 100 5.2.- Métodos de generación de mallados no estructurados...................................... 102 5.3.- Mallados “Multiblock” ..................................................................................... 103 5.4.- Mallados ajustados a los contornos (“Body Fitted Coordinates” o BFC) ........ 1056.- Bibliografía .................................................................... 109Anexo I.- El metodo Von Neumann para analisis deestabilidad ........................................................................... 113Anexo II.- Ecuaciones en sistemas generalizados............... 123Anexo III.- Glosario de términos empleados en técnicasnuméricas ............................................................................ 129 5. Técnicas Numéricas en Mecánica de Fluidos5Prólogo.Las técnicas numéricas en Ingeniería han experimentado un gran desarrollo en lasúltimas décadas. De esta tendencia no se ha apartado una rama tan característica de laIngeniería Mecánica como es la Mecánica de Fluidos. En este trabajo se intenta recopilar lastécnicas más utilizadas y los tratamientos que se pueden hacer de las ecuaciones deconstitución con el fin de obtener una solución numérica de las mismas.Esta publicación constituye el resultado de varios años de impartición de la asignaturade doctorado “Dinámica de Fluidos Computacional” que los profesores del Área de Mecánicade Fluidos de la Universidad de Oviedo han realizado en distintos Departamentos de lamisma. Los autores han tratado de recoger la tradición y conocimientos que comparten con elresto de compañeros del Área, a los que agradecen sus enriquecedores comentarios. 6. 6 Técnicas Numéricas en Mecánica de Fluidos 7. Técnicas Numéricas en Mecánica de Fluidos 71.- Introducción. La dinámica de fluidos computacional (o CFD, acrónimo de las palabras inglesas“Computational Fluid Dynamics”) consiste en el análisis del movimiento de los fluidosmediante simulaciones con ordenadores. Su objetivo es la búsqueda de una soluciónaproximada de las ecuaciones que gobiernan el movimiento de los fluidos, discretizando odividiendo el dominio de cálculo en pequeños elementos y resolviendo allí dichas ecuaciones.Los métodos numéricos aplicados a la Mecánica de Fluidos resultan una herramientamuy útil para el diseño y análisis de las distintas situaciones prácticas en las que se utilizanlos fluidos.1.1.- Perspectiva histórica.Hasta finales de los años 60 los ordenadores no alcanzaron velocidades de cálculosuficientes como para resolver casos sencillos. Hasta entonces, las técnicas experimentalesconstituían prácticamente la única herramienta de análisis y diseño de cualquier problema deMecánica de Fluidos. En la actualidad, los ensayos experimentales siguen siendo necesarios para lacomprobación de las prestaciones de diseños complejos, pero los continuos avances en lacapacidad de los ordenadores y en los algoritmos permiten una reducción importante en elnúmero de ensayos necesarios. Así, por ejemplo, el diseño típico de un modelo de ala de avión,exige en la actualidad 3 o 4 ensayos en túnel aerodinámico, en lugar de los 10 o 15 que erannecesarios anteriormente. En realidad, la palabra que mejor definiría hoy la relación entreambas herramientas podría ser la de complementariedad (Strazisar, 1994, Lakshminarayana,1991).A lo largo de los últimos veinte años, las técnicas de CFD han evolucionado, mejorandolos programas comerciales e introduciéndose en las distintas áreas de la ingeniería hasta hacerseun hueco dentro de las necesidades reales de la industria. Dichos programas se vienen usandode manera creciente paralelalelamente a la mejora de los sistemas informáticos que les sirven desoporte.Desde su inicial concepción, orientada a la industria aeroespacial, la técnicas numéricasse han extendido a un número creciente de aplicaciones dentro de un amplio espectro deindustrias, desde las más clásicas como la automovilística o electrónica, hasta las nuevasaplicaciones en industrias alimentarias y biomédicas. Sin embargo, aún no se emplean lastécnicas numéricas como auténticas herramientas de diseño. En realidad, se obtendrán losmayores beneficios cuando se llegue a un nivel de utilización en el día a día del diseñoindustrial, es decir, cuando su difusión sea tal que puedan llegar a ser utilizadas por ingenierossin demasiada especialización en estas técnicas y no estar restringida su utilización a losexpertos en la materia. 8. 8Técnicas Numéricas en Mecánica de Fluidos1.2.- Aplicaciones habituales de las técnicas numéricas.Es innegable que la industria aeroespacial fue pionera en el trabajo con técnicas CFD yque todavía hoy se encuentra entre las vanguardias en la explotación de estas técnicas, pero cadadía resultan más comunes las aplicaciones en “procesos industriales”, donde los flujos son abaja velocidad y muchas veces prácticamente incompresibles. Esta denominación de procesoindustrial se usa a menudo en la bibliografía técnica (Hirsch, 1995) para distinguir dichasaplicaciones de las aeronáuticas o de flujos a altas velocidades, con grandes efectos decompresibilidad. Entre las aplicaciones más importantes en que se emplean las técnicasnuméricas, se tienen:a) Industria automovilística. Las aplicaciones típicas son el estudio de la aerodinámicade vehículos, la climatización del habitáculo interior, el enfriamiento del bloque motor, el flujoen válvulas de distribución, el diseño de filtros y elementos de control y las investigacionessobre la descarga de combustible en depósitos. b) Industria electrónica. Los problemas más estudiados son el flujo y la distribución detemperaturas en las carcasas electrónicas, el enfriamiento de distintos componentes, el flujo deaire en las unidades de discos, los procesos de construcción de chips usando la técnica dedeposición química del vapor (CVD) y algunos problemas indirectos, como la ergonomía degrandes salas.c) Industrias de proceso y químicas. Problemas habituales resueltos con técnicas CFDson el flujo de plásticos, los estudios en conducción de lodos, el flujo del vidrio fundido, losflujos de tintes, la deposición de vapores químicos, el llenado de moldes, estudios en procesosde combustión y los flujos reactivos complejos (con intercambio de calor, masa y reaccionesquímicas). d) Industria de conformados metálicos. Las aplicaciones más comunes en esta industriason los procesos de fundición continua, las fundiciones abiertas, la extrusión de metales y losprocesos de solidificación (construcción de hélices de barco, por ejemplo).e) Industria nuclear. Algunos estudios relacionados con el flujo en conductos desustancias originadas en los procesos de reacción nuclear, el enfriamiento del reactor, estudiosrelacionados con el intercambiador de calor, el flujo en el interior del reactor, elalmacenamiento de residuos nucleares, el diseño de torres de enfriamiento y las investigacionessobre chorros térmicos.f) Industria de recubrimientos de película fina. Entre otros, los problemas estudiados pormedio de técnicas CFD han sido el recubrimiento de cintas magnéticas, el recubrimiento depelículas de fotografía o de sonido, el recubrimiento de adhesivos, multitud de aplicaciones enla industria papelera y los recubrimientos de fibra óptica. g) Industria biomédica y farmacéutica. Entre otras aplicaciones, destacan el flujo de lasangre en la venas y arterias, el flujo a través de distintas prótesis, el flujo en el interior delcorazón, los distintos estudios en fenómenos de centrifugación y el diseño de sistemas deinyección intravenosa. 9. Técnicas Numéricas en Mecánica de Fluidos9 h) Industria alimentaria. Destacan los diseños de procesos de pasteurización, los estudiosen equipos de procesado de alimentos con estructura toroidal, la extrusión de fluidos y loshornos de convección. i) Industria aeroespacial. Las aplicaciones habitualmente estudiadas son los efectos de lamicrogravedad, la ventilación de habitáculos, el diseño de vehículos espaciales, los flujos decombustible en conductos y tanques, estudios varios en motores de propulsión.j) Industria aeronáutica y naval. Estudios en perfiles aerodinámicos, diseño de trenes deaterrizaje, estudios en hélices marinas y el diseño de carenas de barcos.k) Otras aplicaciones. Destacan los estudios en oceanografía, predicciones en hidrología(planificación de embalses, regímenes de precipitaciones, entre otros), los flujo en conductos(calefacción, flujos internos en edificaciones, ingeniería de complejos urbanos), los flujosmedioambientales, la meteorología, los estudios de flujos alrededor de edificios, puentes y otrasestructuras exteriores, las investigaciones relacionadas con la propagación de humos, losestudios sobre el enfriamiento y crecimiento del vidrio, el flujo en máquinas de desplazamientopositivo, las investigaciones en flujos con varias fases (sprays), el estudio de los MEMS (MicroElectro-Mechanical Systems) y las aplicaciones en turbomáquinas.Por lo tanto, la panorámica es realmente amplia y susceptible de crecimiento en elfuturo.Las ventajas que proporciona el análisis con técnicas CFD se pueden resumir en lassiguientes:- Reducción sustancial de tiempos y de costes en los nuevos diseños.- Posibilidad de analizar sistemas y condiciones muy difíciles de simularexperimentalmente: velocidades supersónicas, temperaturas extremas y elementosen movimiento relativo.- Capacidad de estudiar sistemas bajo condiciones peligrosas o más allá de suscondiciones límite de funcionamiento, por ejemplo, accidentes con sustanciastóxicas.- Nivel de detalle prácticamente ilimitado. Los métodos experimentales son tantomás caros cuanto mayor es el número de puntos de medida, mientras que losprogramas CFD pueden generar un gran volumen de resultados sin coste añadido,y con posibilidad hacer estudios paramétricos.- Un valor añadido es poder poner en el producto la etiqueta de “Diseñado con ayudadel ordenador”, y la facilidad para generar gráficos fácilmente interpretables, queestimulan la “compra” del producto (lo cual, por otro lado, constituye un riesgo).Aunque se razone que se abaratan los costes respecto a las técnicas experimentales, lastécnicas CFD no son gratuitas. En primer lugar, son necesarias máquinas de gran capacidadde cálculo (los usuarios de técnicas CFD trabajan habitualmente con los ordenadores máspotentes que existen), y un software con precio todavía no accesible al gran público. Ensegundo lugar, se necesita personal cualificado que sea capaz de dominar los programas y,sobre todo, analizar adecuadamente los resultados. Los desarrollos en el campo de las técnicas numéricas dedicadas al estudio de flujos seestán acercando cada vez más a los de otras herramientas de CAE como las de análisis de 10. 10 Técnicas Numéricas en Mecánica de Fluidosesfuerzos en sólidos y estructuras. Su retraso se debe a la gran complejidad de las ecuacionesy, sobre todo, a la dificultad de modelizar adecuadamente ciertos fenómenos como laturbulencia, los flujos multifásicos y la combustión.Uno de los mayores inconvenientes de las técnicas CFD consiste en que no siempre esposible llegar a obtener resultados suficientemente precisos y siempre está presente laposibilidad de cometer graves errores en cuestiones básicas. Esto proviene de: - Simplificación del fenómeno a estudiar para que el hardware y software sean capaz de abordarlo. El resultado será tanto más preciso cuanto más adecuadas hayan sido las hipótesis y las simplificaciones realizadas. - La existencia de insuficientes e incompletos modelos para la simular el efecto de la turbulencia, los flujos multifásicos o la combustión, entre otros. - La tendencia humana de creerse todo lo que se ha obtenido utilizando un ordenador, sobre todo cuando se presentan los resultados en forma atractiva.1.3.- Niveles de aproximación empleados en las técnicas numéricas.El desarrollo de las técnicas numéricas y su aplicación a cualquier ciencia o tecnologíahan dado lugar al desarrollo y a la concienciación generalizada de uno de los conceptos básicosen ingeniería como es el de grado de aproximación. Esta idea es bastante clara si se consideraque lo que se pretende con cualquier técnica numérica es conocer las variables físicas a partir dela resolución numérica de una serie de ecuaciones que gobiernan el fenómeno.Se han de definir y establecer las distintas aproximaciones que introducen los métodosnuméricos. En lo referente a la Mecánica de Fluidos, la primera aproximación que aparece es elplanteamiento del modelo físico-matemático que defina el comportamiento real de undeterminado flujo. Dicho modelo matemático está habitualmente basado en la hipótesis delcontinuo, válida para la mayor parte de problemas industriales, pero que tiene sus limitacionespara casos extremos de flujos de gases. Una vez hecha esta salvedad, aplicando las leyes básicasde la física clásica se puede establecer una serie de ecuaciones diferenciales que gobiernan elcomportamiento matemático de toda partícula fluida. La resolución exacta de dichas ecuacionesserviría para determinar completamente cualquier movimiento en el seno de un fluido. Se puededecir que un modelo matemático se define únicamente tras haber considerado el nivel deaproximación a la realidad requerido a la hora de obtener la exactitud deseada en el cálculo deuna serie de variables dependientes. Desafortunadamente, debido a la complejidad de lasecuaciones diferenciales que aparecen, a la complejidad geométrica de los flujos, y a lacomplejidad de las condiciones de contorno o iniciales, no resulta posible obtener solucionesanalíticas de dichas ecuaciones de gobierno.Establecidas las ecuaciones de gobierno resulta imprescindible introducir una segundaaproximación al problema. La forma clásica de abordarlo sería construir un modelo a escalareducida del flujo en cuestión y analizarlo experimentalmente en el laboratorio. Laaproximación numérica implica introducir algunas hipótesis simplificativas que aproximen lomás posible los resultados finales a los que se obtendrían si se pudiera calcular la soluciónexacta. Dichas hipótesis se dirigen habitualmente hacia la simplificación tanto de la geometría aestudiar como de las ecuaciones a resolver. Obviamente, al no disponerse de la soluciónanalítica exacta resulta bastante complicado establecer de antemano qué hipótesis sirven y 11. Técnicas Numéricas en Mecánica de Fluidos 11cuales son descartables y, por tanto, en cualquier simulación aplicada a la Mecánica de Fluidos,es preciso dedicar mucho esfuerzo al análisis de los resultados obtenidos antes de aceptarloscomo válidos.Una vez definidas las ecuaciones diferenciales simplificadas, aparece otro problemarelacionado con el posible tratamiento que se pueda hacer de dichas ecuaciones usando técnicascomputacionales. Por medio de los ordenadores resulta muy fácil resolver una ecuación osistema de ecuaciones algebraico, sin embargo, las ecuaciones que estudian el movimiento delos fluidos son ecuaciones diferenciales no lineales. Resulta obligatorio realizar latransformación de las ecuaciones de forma que puedan ser resueltas por un ordenador. El pasode las ecuaciones diferenciales a sus equivalentes lineales constituye otro nivel de aproximacióny normalmente recibe el nombre de discretización de las ecuaciones.En cuanto a la geometría a estudiar, se debe señalar que la aproximación a la que debesometerse no sólo es de orden descriptivo respecto a su contorno sino que además ha deestablecerse la definición del espacio ocupado por el fluido. En este sentido, resultaimprescindible referir los puntos a un determinado sistema de coordenadas en los que sepretenderá resolver las ecuaciones para obtener soluciones de las variables deseadas. Aunque elcampo fluido sea un continuo, no se puede pretender resolver las ecuaciones en todos los puntosde un determinado volumen, porque entonces se tendría un número enorme de ecuaciones aresolver. Por tanto, hay que elegir cierto conjunto de puntos en los que se resolverán lasmencionadas ecuaciones y que serán los puntos dónde finalmente se conocerán los valores delas variables fluidas. La definición de estos puntos es lo que se denomina habitualmentediscretización espacial del dominio (también se habla de generación del mallado). El procesodescrito no deja de ser otra aproximación que se introduce en el cálculo y que define el nivel deaproximación espacial. En el caso de tener ecuaciones que dependan de la variable tiempo (flujo noestacionario) es esencial la definición de un nivel de aproximación temporal. No es posibletampoco estudiar la evolución de las variables en el tiempo de forma continua. El nivel indicarála forma de modelizar la evolución real introduciendo lo que se denomina discretizacióntemporal del sistema de ecuaciones. A partir de la solución calculada se podrá realizar unpromediado temporal oportuno para estudiar ciertas características medias del flujo quedependan de la evolución de las variables con el tiempo.Finalmente, se pueden manipular las ecuaciones eliminando ciertos términos cuyainfluencia en un determinado problema se considere despreciable. La conclusión de que algúntérmino no afecta a la solución de una determinado flujo se debe alcanzar tras analizardetenidamente la sensibilidad del problema ante valores dispares de dicho término.Normalmente dicho estudio se hace tras adimensionalizar convenientemente las ecuaciones yrealizar el correspondiente análisis de semejanza (técnicas asintóticas). Esta cuestión es deimportancia capital en la Mecánica de Fluidos y está en el origen de cualquier estudioexperimental. Desde el punto de vista numérico, la eliminación de algún término en lasecuaciones introduce lo que se denomina nivel de aproximación dinámico de las ecuacionesconsideradas. Resumiendo, desde el modelo matemático (ecuaciones diferenciales no lineales) queaproxima la realidad física en un medio continuo se llega a un número finito de ecuacionesalgebraicas que eliminan algún término de las ecuaciones de partida y que aproximan laevolución temporal real que, tras resolver con técnicas apropiadas, proporcionan una 12. 12Técnicas Numéricas en Mecánica de Fluidosaproximación al valor de las variables incógnita en los puntos elegidos como discretizaciónespacial del dominio de cálculo. En definitiva, se establece un nivel de aproximación numéricolímite por debajo del cuál será imposible acercarse al valor real de las variables en los puntoselegidos. Sin embargo, desde un punto de vista ingenieril, el proceso descrito es perfectamenteválido y ha significado a lo largo de la evolución de las técnicas numéricas, la posibilidad demejorar diseños y ahorrar mucho esfuerzo que de otra manera supondría trabas insalvables a laevolución de muchos sectores industriales. En la figura 1 se muestra gráficamente lapanorámica explicada en este apartado.Figura 1.- Las técnicas numéricas en la Mecánica de Fluidos. 13. Técnicas Numéricas en Mecánica de Fluidos132.- Modelos físico-matemáticos.Existen varias posibles simplificaciones en cuanto a la definición del modelomatemático que describe el movimiento de las partículas de un fluido, de gran interés por darlugar a soluciones válidas en distintos problemas, que han sido ampliamente utilizadas enmuchas aplicaciones numéricas. Se enumeran y describen brevemente los modelos másimportantes. Normalmente cada una de estas estrategias ha sido desarrollada para un tipoparticular de flujo.2.1.- Modelos de flujo potencial.Describen el comportamiento de flujos irrotacionales de fluidos ideales, desarrolladosen los albores de las técnicas numéricas. La teoría básica para el cálculo consiste en partir dela definición del concepto de potencial de velocidades. Constituye una simplificaciónadicional muy empleada para el cálculo de flujos estacionarios. Conceptualmente es de graninterés, pero está cayendo en desuso.2.2.- Ecuaciones para flujo ideal.Cuando el número de Reynolds es suficientemente elevado, lo que ocurre en muchasde las aplicaciones prácticas de la Mecánica de Fluidos, despreciar los efectos viscosos y deconducción de las ecuaciones resulta una aproximación bastante cómoda pues elimina lostérminos difusivos de segundo orden en las ecuaciones diferenciales y hace que las ecuacionesde gobierno pasen a ser de primer orden, con todo lo que lleva asociado en cuanto asimplificación de cálculos. Con las hipótesis de despreciar los efectos viscosos y la transferencia de calor porconducción, es decir, si se considera al fluido como ideal, se obtienen las ecuaciones de Euler.Los modelos numéricos que resolvían las ecuaciones de Euler eran hasta hace poco los únicosexistentes. Hoy en día este tipo de modelos constituye el punto de partida para el desarrollo demodelos más completos (Arts, 1989).Las ecuaciones de Euler fueron desarrolladas por este famoso matemático suizo haciael año 1670. Adoptan las expresiones siguientes: Continuidad:dρ + ρ∇·u = 0(1)dt 14. 14Técnicas Numéricas en Mecánica de Fluidos Cantidad de movimiento:∂uρ+ ρ(u·∇)u = −∇p + ρ f e (2)∂tdonde f e es la resultante de las fuerzas volumétricas externas que afecta a cada partícula.Habitualmente sólo se considera la fuerza debida a la gravedad, pero puede tener otroscomponentes: fuerzas electromagnéticas, fuerzas de Coriolis, fuerzas centrífugas, etc.2.3.- Modelo para flujo ideal, estacionario y rotacional. Es un tipo de modelos muy similar a los del fluido ideal. Consiste en reducir elnúmero de variables que intervienen en los cálculos introduciendo la vorticidad en lasecuaciones de cantidad de movimiento y de la energía. Normalmente se parte de ladenominada representación de Clebsch para la velocidad en función de la vorticidad. No seconsideran aquí ni las pérdidas por viscosidad en la capa límite ni los efectos de laturbulencia.2.4.- Solución directa de las ecuaciones de Navier-Stokes.Las ecuaciones de Navier-Stokes constituyen una modelización correcta del flujo deun fluido Newtoniano, incluyendo todos los efectos viscosos y térmicos. Adecuadamenteresueltas incluyen los efectos de la turbulencia y de la capa límite. Pero esta resolución directa(ver Hirsch, 1988) requiere de una discretización espacial y temporal tan fina que estáclaramente fuera del alcance de cualquier aplicación industrial. Se ha estimado (Moin et al.,1997) en varios cientos de años de CPU del ordenador más potente existente en aquel año elcálculo de un segundo de vuelo de un avión comercial.La resolución numérica de las ecuaciones de Navier-Stokes sí que es posible si seutilizan modelos adecuados para simular el efecto de la turbulencia y de la capa límite endiscretizaciones no tan detalladas. Un poco más adelante se hablará de ellos.2.5.- Modelo parabólico de las ecuaciones de Navier-Stokes (PNS). Este tipo de modelos ha sido desarrollado para el cálculo de flujos supersónicos ehipersónicos, donde la captura de las ondas de choque, gradientes de presión, esfuerzosviscosos superficiales y transferencia de calor son los objetivos más importantes paracualquier diseño. Las ecuaciones de gobierno parabólicas se obtiene a partir de las de Navier-Stokesconsiderando las siguientes hipótesis: - Flujo estacionario. 15. Técnicas Numéricas en Mecánica de Fluidos 15 - Los gradientes de esfuerzos viscosos son despreciables en la dirección de las líneas de corriente. - Los gradientes de presión en la dirección de las líneas de corriente se aproximan por su valor en zonas de capa límite cercanas. Se ha investigado mucho en las ecuaciones parabólicas de Navier-Stokes (ParabolizedNavier-Stokes) y se han desarrollado muchos algoritmos en los últimos años (Hirsch, 1988 yHoffman, 1989), aunque su aplicación está muy limitada al sector aeroespacial.2.6.- Modelo de flujo incompresible. Un flujo se denomina incompresible cuando la densidad del fluido en cada instantepermanece independiente de las variaciones de presión. La importancia de los flujosincompresibles es indudable y algunos autores, como Batchelor (1967), llegan a afirmar que losproblemas relacionados con este tipo de flujos constituyen la aplicación más importante ycompleja de resolver de la Mecánica de Fluidos.Cuando el flujo es además isotermo, las ecuaciones de gobierno se simplificannotablemente y la solución para las distintas variables se hace independiente de latemperatura. El sistema de ecuaciones requerido queda reducido a la ecuación de continuidady la de cantidad de movimiento, que expresadas adimensionalmente y con la única presenciade la gravedad como fuerza volumétrica, adoptan la forma: Continuidad:∇·u = 0 (3) Cantidad de movimiento:du∇p 1 21 =- +∇ u+g(4)dtρ ReFr Contrariamente a lo que pudiera pensarse, la hipótesis de incompresibilidad complicabastante la resolución de las ecuaciones. No sólo la densidad sino también los distintoscoeficientes de transporte del fluido son independientes de la presión y de la temperatura. De estaforma, las ecuaciones de continuidad y de cantidad de movimiento son independientes de laecuación de la energía, que no es necesario resolver para obtener los campos de velocidades ypresión. Pese a la ventaja que esto parece implicar, en la práctica, las dos ecuaciones a resolver sevuelven “rígidas” por la ausencia de derivada temporal en la ecuación de continuidad, y susolución resulta más dificultosa al no ser posible una iteración tomando ambas como punto departida.2.7.- Modelos para simulación de la turbulencia. El número de Reynolds de un flujo da una medida de la importancia relativa de lasfuerzas de inercia, asociadas con los efectos convectivos, y las fuerzas viscosas. En 16. 16 Técnicas Numéricas en Mecánica de Fluidosexperimentos con fluidos se observa que para valores inferiores a un número de Reynoldsdenominado crítico, el flujo es intrínsecamente estable y las capas de fluido adyacentes sedeslizan unas sobre otras de forma ordenada. El régimen del flujo se denomina laminar.Si el flujo tiene un valor del número de Reynolds por encima del denominado crítico,se manifiestan en éste unas perturbaciones que dan lugar a un cambio radical del carácter delflujo. El movimiento se vuelve intrínsecamente no estacionario, incluso con condiciones decontorno constantes. Este régimen se denomina flujo turbulento.La turbulencia se define como el estado de movimiento de un fluido en el que lasdistintas variables relevantes (presión, velocidad, etc.) fluctúan de una forma desordenada. Setrata de un estado no estacionario desde el punto de vista macroscópico en el que las distintasvariables adoptan valores dependientes tanto de la posición como del tiempo y estos valoresvarían de una forma aleatoria y desordenada. La descripción del movimiento de las partículas fluidas debido al efecto de laturbulencia resulta altamente complejo y constituye un problema aún sin solución desde elpunto de vista de los métodos numéricos. Se han propuesto varias formas de resolver elproblema utilizando distintas aproximaciones. A continuación se exponen los métodosconocidos como simulación directa, simulación de los grandes vórtices y promediadotemporal de las ecuaciones de Navier-Stokes. - Simulación directa de las ecuaciones (“Direct Simulation”, DS).Este método (cuyas iniciales provienen de la denominación inglesa Direct Simulation)consiste, en realidad, en no utilizar ningún modelo para la turbulencia, sino realizardiscretizaciones temporales y espaciales que sean capaces de simular el flujo en undeterminado problema.La resolución directa de las ecuaciones de Navier-Stokes resulta hoy en día abordablesólo para un número muy limitado de problemas simples de interés académico. Los grandescentros dedicados a la Mecánica de Fluidos disponen de líneas de investigación con estaorientación, pero tanto las limitaciones en memoria de almacenamiento de las variables, comoen el tiempo de cálculo hacen de momento impensable la solución generalizada de problemasprácticos usando este tipo de técnicas. Según Vandromme (1989), la primera solución de estetipo se realizó en 1981 en la Universidad de Stanford. - Simulación de grandes vórtices (“Large Eddy Simulation”, LES). Este tipo de técnicas numéricas reducen la complejidad de las ecuaciones de gobiernoconsiderando sólo parte de los efectos turbulentos del flujo. Se estudia el intercambioenergético entre las denominadas “fluctuaciones de gran escala” y se simula el efecto de laspequeñas escalas de la turbulencia.Se trata de un tipo de modelo intermedio entre la simulación directa y el promediadotemporal de las ecuaciones de Navier-Stokes, que extiende el promedio temporal a la capturade ciertos efectos turbulentos básicos de forma numérica. En los modelos de simulación de 17. Técnicas Numéricas en Mecánica de Fluidos17grandes vórtices, las ecuaciones no estacionarias del flujo se resuelven para el flujo principaly para los vórtices grandes mientras que se modela el efecto de los vórtices pequeños. Aunque sin llegar al extremo de la simulación directa, sólo es posible para problemassimplificados y requiere unas capacidades de cálculo muy elevadas. - Modelos que promedian temporalmente las ecuaciones de Navier-Stokes (RANS). Los modelos de promedio de las ecuaciones de Navier-Stokes (Reynolds AveragedNavier-Stokes) han sido muy estudiados y resultan bastante útiles en la mayoría de losproblemas prácticos resueltos mediante técnicas numéricas. El procedimiento de promediar las leyes que describen el movimiento de una partículase introduce en las ecuaciones con el fin de obtener los comportamientos promedio yturbulento (aleatorio) de las distintas variables. El punto de partida es muy sencillo. Se tratade obtener una descomposición de las variables en su valor medio y su valor fluctuante. Porejemplo, para la velocidad, la descomposición sería:u = u + u (5)donde la componente media de la velocidad se obtiene haciendo la integral de la velocidadinstantánea: 1 T T∫0 u (t) = u(t) dt (6)suponiéndose que el periodo de integración (T) es lo suficientemente grande en comparacióncon la escala temporal de la turbulencia, pero lo suficientemente pequeño como para captarcualquier fenómeno no estacionario distinto a la turbulencia. La utilización de este tipo demétodos es bastante adecuado, pues la mayoría de los fenómenos no estacionarios enMecánica de Fluidos tiene lugar a frecuencias con rangos muy alejados del rango defrecuencias de la turbulencia.El proceso de promediado temporal de las ecuaciones diferenciales, da lugar a unostérminos, denominados de tensiones de Reynolds (“Reynolds stresses”), que involucranmedias de los productos de la fluctuaciones de las componentes de la velocidad, cuya relacióncon las componentes medias del flujo es desconocida. Para obtener dicha relación esnecesario introducir un modelo adicional, denominado modelo de turbulencia o de cierre. Lasdistintas posibilidades prácticas en cuanto a modelos de turbulencia son analizadas acontinuación. Habitualmente lo que interesa son los efectos de la turbulencia sobre los valoresmedios de las variables: la velocidad media y la presión media en el caso del flujo en unconducto; en el caso de un avión, las fuerzas medias de resistencia y sustentación; para el casode un motor, los efectos de la turbulencia sobre las relaciones de mezcla entre combustible ycomburente. 18. 18Técnicas Numéricas en Mecánica de FluidosPara conseguir esto, las ecuaciones de Navier-Stokes se promedian sobre las escalasde las fluctuaciones de turbulencia (RANS). Estos métodos dan lugar a un campo de flujopromediado y simulado que es más uniforme que el flujo real, y, por tanto, reducedrásticamente el número de puntos de la discretización espacial y de la temporal necesariopara obtener las variables buscadas. Un modelo de turbulencia es un procedimiento numérico que permite relacionar losvalores medios de las fluctuaciones de las variables con los valores promedio (en lanomenclatura propia de estos métodos, se habla de cerrar el sistema de ecuaciones), de formaque se puedan resolver las ecuaciones de gobierno. Un modelo de turbulencia será útil, dentrode un programa CFD de propósito general, si es exacto, sencillo y económico(computacionalmente hablando). Los modelos de turbulencia más comunes son los siguientes(Wilcox, 1993): Modelos algebraicos: Modelo de cero ecuaciones: modelo de la longitud de mezcla. Modelo Cebeci-Smith-Mosinki. Modelo Baldwin-Lomax. Modelos de viscosidad turbulenta: Modelo de una ecuación: modelo k. Modelo de dos ecuaciones: modelos k-ε, k-ω2 o q-ω, modelo RNG. Modelos de ecuaciones de las tensiones de Reynolds (RSM).Cada uno de ellos tiene sus ventajas e inconvenientes. Aquí se expondrá más endetalle el modelo k-ε, por ser muy extendido.2.8.- Modelos de dos ecuaciones: modelo k-ε.Como modelo de cierre o estrategia numérica para resolver de forma aproximada lasecuaciones de Navier-Stokes, se desarrollan dos ecuaciones de transporte adicionales(similares a las definidas por las ecuaciones 3 y 4), una para la energía cinética turbulenta (k)y otra para la tasa de disipación de energía cinética turbulenta (ε). Estas variables se definensegún las expresiones: k=2(1 2u+ v 2 + w2 )(7)ε = 2 ν e ij e ij(8)donde e ij es la parte fluctuante del tensor de velocidad de deformación. Las ecuaciones de transporte para k y ε se basan en el conocimiento de los procesosque producen los cambios en esas variables y son (Versteeg et al., 1995): 19. Técnicas Numéricas en Mecánica de Fluidos19∂ ( ρk ) µ  + ∇·( ρku ) = ∇· t grad k  + 2µ t E ij E ij − ρε(9)∂t  σk  ∂ (ρ ε ) µ ε ε+ ∇·( ρ ε u ) = ∇· t grad ε  + C1ε 2 µ t E ij E ij − C2ε ρ ε(10) ∂t  σε k kdonde Eij es el tensor de componentes medias de la velocidad de deformación. El significadofísico de las anteriores expresiones se puede resumir en el siguiente balance: Velocidad de  Transporte de k / ε  Transporte de k / ε   Producción   Destrucción += + de k / ε − de k / ε cambio de k / ε   por convección   por difusión    Aparecen varios conceptos cinemáticos relacionados con las “escalas” o longitudestípicas asociadas a los distintos movimientos del flujo (flujo principal medio y flujo oscilanteo turbulento, relacionado con los vórtices). La escala de velocidad ζ es característica de losremolinos y de las propiedades del flujo principal y se define según la expresión:∂ u ζ =C (11)∂ ydonde C es una constante adimensional yla escala de longitud turbulenta (o longitud demezcla), que se define como: k 3/ 2 =(12) εEste método utiliza la velocidad de disipación ε de los remolinos pequeños paradefinir la escala de longitud de los remolinos grandes porque, para altos números deReynolds, la velocidad de extracción de energía del flujo de los remolinos grandes es igual ala velocidad de transferencia de energía a los remolinos pequeños. Si esto no fuese así, laenergía en algunas escalas de la turbulencia podría aumentar o disminuir sin límite, cosa queno ocurre en la práctica con lo que se justifica el uso de la velocidad de disipación ε dentro dela definición de la escala de longitud “ ”. Aplicando la misma aproximación del modelo de la longitud de mezcla se puedeobtener la viscosidad turbulenta como: k2µ t = ρC µ ζ = ρC µ (13) ε Las ecuaciones de transporte de k y ε contienen cinco constantes ajustables Cµ, losnúmeros de Prandt (σk y σε), C1ε y C2ε, aunque se suelen emplear valores fijos para unaamplia gama de flujos turbulentos. Los números de Prandtl (números adimensionales quemuestran el peso relativo de los términos viscosos frente a los términos de transmisión decalor por conducción) relacionan las difusividades de k y ε con la viscosidad turbulenta (µt). 20. 20Técnicas Numéricas en Mecánica de Fluidos Las principales ventajas e inconvenientes del modelo, tal y como pueden encontrarseen la bibliografía consultada (Lakshminarayana, 1991 y Versteeg et al., 1995), son lassiguientes:a) Ventajas: - Sólo se necesita fijar las condiciones iniciales y de contorno. - Resultados satisfactorios para una gran cantidad de flujos. - Es el modelo turbulento más ampliamente utilizado en la mayoría de flujos en aplicaciones industriales. - Se dispone de leyes de pared desarrolladas como condiciones de contorno para este tipo de modelos.b) Inconvenientes: - Implementación más compleja que los modelos algebraicos debido a la introducción de dos ecuaciones diferenciales adicionales. - Pobres resultados en casos como: flujos no confinados, flujos con grandes gradientes longitudinales, flujos turbulentos completamente desarrollados en conductos no circulares.2.9.- Modelos matemáticos para la Capa Límite.La capa límite es la zona del campo fluido próxima a un contorno sólido en la que semanifiestan especialmente los efectos viscosos. Debido a la viscosidad y a la condición de nodeslizamiento, cerca de cualquier contorno sólido aparece un gradiente de velocidades en ladirección normal a dicho contorno. Este gradiente de velocidades condiciona el intercambioenergético entre las distintas partículas de fluido con velocidades diferentes, originandovorticidad y turbulencia.El problema básico para la modelización numérica del intercambio energético en lacapa límite sobre cualquier frontera sólida consiste en la definición correcta de las velocidadesde las partículas en una zona muy próxima a dicha frontera. Esto implica una densidad demallado muy elevada, necesaria para capturar los distintos fenómenos que se producen dentrode la capa límite. Esta dificultad se ha abordado usando varias aproximaciones, que se pueden englobaren cuatro grupos: modelos de distribución de las pérdidas, modelos de capa de cortadura,modelos de capa límite y leyes de pared, que son brevemente explicados a continuación. - Modelos de distribución de las pérdidas (“Distributed Loss Models”). Este tipo de modelos constituye una aproximación muy usada en flujos internos (elfluido está confinado en un canal de paso limitado por paredes sólidas). La hipótesis básicaconsiste en suponer que el efecto de las tensiones cortantes debidas a la viscosidad es 21. Técnicas Numéricas en Mecánica de Fluidos 21equivalente a una fuerza de rozamiento distribuida a lo largo del canal de paso y definida porvalores semi-empíricos conocidos del problema a resolver. Aunque con este tipo de modelos se puede predecir el flujo en gran parte de lageometría, es claro que se pierde la definición en zonas cercanas a las superficies sólidas. Aveces esta falta de precisión en la definición del flujo no es tolerable y se requiere superponeralgún modelo de capa límite complementario.Los modelos de distribución de pérdidas fueron muy populares en los inicios de lastécnicas CFD cuando la potencia de cálculo hacía difícil de llevar a la práctica cualquier otrotipo de modelo (Bosman, 1976). - Modelos de capa límite (“Boundary Layer Approximations”). Derivado de los estudios de Prandtl sobre la estructura del flujo para elevados valoresdel número de Reynolds. Bajo estas condiciones, el campo de velocidades en un fluido sepuede separar en dos zonas, una de flujo no viscoso alejada de los contornos sólidos y otradominada por los efectos de la viscosidad, denominada capa límite.Las ecuaciones de este tipo de modelos se pueden derivar de las del modelo de la capade cortadura simplificándolas aún más mediante la hipótesis del valor despreciable de lavelocidad en la dirección normal a la superficie considerada en comparación con la velocidaden la dirección de las líneas de corriente.También existen muchas aplicaciones prácticas de este tipo de modelos (Launder etal., 1972) - Modelos de la capa de cortadura (“Thin Shear Layer, TSL”). Son métodos apropiados para flujos con elevados números de Reynolds en los que laszonas de influencia viscosa, estelas o capas de cortadura ocupan una extensión muy reducidadentro de la geometría del problema estudiado. Fuera de estas zonas, resulta suficiente conconsiderar el modelo de fluido ideal. Para este tipo de modelos se requiere una discretización espacial muy densa en laszonas en las que se espera influencia de los términos viscosos. En realidad, se trata de uncálculo ligeramente más avanzado que el correspondiente al modelo de capa límite, porque eneste caso la geometría de la capa límite es resultado del cálculo y no se introducen hipótesisadicionales. Este tipo de modelos ha sido aplicado a multitud de problemas relacionados conaplicaciones aerodinámicas (Hirsch, 1988). - Leyes de pared. Una posibilidad distinta a los modelos mencionados consiste en incluir en los cálculosalguna aproximación para la distribución de velocidades esperada. Con tal fin, se pueden 22. 22 Técnicas Numéricas en Mecánica de Fluidosutilizar las distribuciones de velocidad obtenidas experimentalmente, pero la práctica habitualconsiste en utilizar los datos de distribuciones teóricas. En el contexto de los métodosnuméricos, las funciones o “leyes de pared” constituyen un conjunto de fórmulas semi-empíricas que relacionan los valores de las distintas variables en las zonas próximas a loscontornos sólidos y sobre dichos contornos. Normalmente incluyen tanto las relaciones paralas variables medias y fórmulas para el tratamiento de la turbulencia en zonas próximas a loscontornos sólidos. La definición de las distintas fórmulas, con rangos de aplicabilidad variables,provienen de los estudios sobre capa límite y parten de la definición de las variablesadimensionales características de dichos estudios. Suelen distinguirse dos zonas que dan lugara la utilización de las denominadas leyes para capas internas y leyes para capas externas.En la bibliografía existen multitud de modelos basados en alguna hipótesis sobre ladistribución de velocidades dentro de la capa límite (Launder et al., 1972 o White, 1991). 23. Técnicas Numéricas en Mecánica de Fluidos 233.- Aspectos matemáticos del procedimientode cálculo. Los modelos matemáticos de la mayoría de fenómenos físicos pueden ser expresadoscomo un sistema de ecuaciones en derivadas parciales cuasi-lineales de primer o segundoorden. Las propiedades matemáticas de estos sistemas de ecuaciones son un reflejo de laspropiedades físicas de los flujos. Las ecuaciones de flujo representan un balance entre fenómenos de convección ydifusión además de la inclusión de términos fuente. En los flujos difusivos aparecen términoscon derivadas de segundo orden como consecuencia de la Ley de Fick generalizada. Por elcontrario en los flujos convectivos aparecen derivadas de primer orden que expresanpropiedades de transporte o arrastre. Cada una de estas contribuciones tiene influencia en lanaturaleza matemática de las ecuaciones. En concreto se distingue entre los siguientes tipos desistemas de ecuaciones en derivadas parciales: • Elípticos. • Parabólicos. • Hiperbólicos. Antes de dar una descripción matemática rigurosa de cada tipo, se considera comoejemplo la componente sobre el eje x de la ecuación de Navier-Stokes para un flujo laminar eincompresible en coordenadas cartesianas: ∂u ∂p+ ρ ( u·∇ ) u = -+ µ ∆u(14) ∂t ∂x Si se adimensionaliza la ecuación 14 utilizando unas magnitudes de referencia L, T yV se obtiene: VT ∂ u∂p 1+ ( u ·∇ ) u = −+ ∆u (15)L ∂t ∂ x Re donde:ρVL VL • Re == es el número de Reynolds y µ νVT • el número de Strouhal. L Se analizan a continuación dos casos particulares: 24. 24 Técnicas Numéricas en Mecánica de Fluidos Para Números de Reynolds muy bajos (un flujo muy viscoso), el término convectivo puede ser despreciado, quedando por tanto:V2 T ∂ u∂p−+ ∆ u = Re(16) ν ∂t ∂x Analizando la ecuación anterior puede verse que, teniendo en cuenta la naturaleza del flujo (estacionario o no estacionario) aparecen dos tipos de ecuación:∂u • Flujo Estacionario ⇒ = 0.∂t∂p∆u = Re(Comportamiento elíptico) (17)∂x ∂u • Flujo No Estacionario ⇒≠ 0. ∂t V2 T ∂ u∂p∆u =+ Re(Comportamiento parabólico)(18)ν ∂t ∂x Para Números de Reynolds muy altos y fuera de la capa límite, los efectos viscosos apenas tienen influencia en el flujo. La ecuación queda entonces reducida a la ecuación de Euler (Ecuación de primer orden).∂u ∂p + ρ ( u·∇ ) u = -(Comportamiento hiperbólico) (19)∂t ∂x La distinción entre unas ecuaciones y otras tiene una gran importancia porque losmétodos numéricos de discretización y resolución de problemas son diferentes en cada caso.Los fenómenos de difusión actúan en todo el espacio independientemente de la direcciónpredominante del flujo (comportamiento elíptico). Por el contrario, los fenómenos deconvección actúan en la dirección de propagación, en regiones concretas del espacio(comportamiento hiperbólico)Entre ambas posibilidades, una ecuación con comportamiento parabólico representa unasituación intermedia que se puede interpretar como un proceso de difusión en todas lasdirecciones pero amortiguado en el tiempo.En las ecuaciones de Navier-Stokes no estacionarias, la ecuación de continuidad eshiperbólica mientras que las ecuaciones de la cantidad de movimiento y de la energía sonparabólicas en espacio y tiempo. Se habla de un comportamiento parabólico-hiperbólico. Sinembargo, en las ecuaciones de Navier-Stokes en su versión estacionaria, las ecuaciones de lacantidad de movimiento y de la energía tienen un comportamiento elíptico por lo que setienen propiedades elíptico-hiperbólicas. 25. Técnicas Numéricas en Mecánica de Fluidos 253.1.- Condiciones iniciales y de contorno. Se dice que un problema cuya resolución depende de una ecuación diferencial enderivadas parciales está “bien planteado” si la solución depende de manera continua de lascondiciones iniciales y de contorno. Se consideran dos tipos de problemas: • Problema de Condición Inicial. Se conoce la solución en t = 0 y se busca la evolución de dicha solución en el tiempo. • Problema de Condición de Contorno. Se fijan las condiciones en los contornos del dominio y se busca la solución en su interior. En este caso existen tres variantes posibles: condición de Dirichlet, condición de Neumann y condición de Robin, que se detallarán en el apartado 4.En este apartado se analiza cómo se transmite la información relativa a la solución enlas distintas regiones de flujo. Se detalla el fenómeno de propagación de la solución para cadauno de los tipos de clasificación de las EDPS. • Sistemas Hiperbólicos. En la figura 2 se considera un problema hiperbólico bidimensional: Γ constituye uncontorno del dominio y la solución U sobre el trozo de contorno AB se propaga en eldominio del flujo siguiendo las superficies características S A , S B que surgen en AB .Existe una región, llamada “zona de Dependencia”, que delimitan las superficiescaracterísticas que surgen de A y de B junto con el trozo de contorno AB y que determinanla solución en el punto P . Asímismo existe una región, llamada “zona de Influencia”, en lacual la solución resulta influenciada por la solución en el punto P .Figura 2.- Zonas de influencia y de dependencia de un problema hiperbólico. 26. 26Técnicas Numéricas en Mecánica de Fluidos • Sistemas Parabólicos. En la figura 3 se muestra lo que ocurre en un problema parabólico bidimensional.Figura 3.- Zonas de influencia y de dependencia de un problema parabólico.En este tipo de problemas ambas superficies características no tienen más que unalínea en común ( S A ). La “zona de Dependencia” del punto P se extiende a un lado de lalínea S A y la “zona de Influencia” ocupa el semiespacio restante. • Sistemas Elípticos. En la figura 4 se considera un problema elíptico bidimensional. En este caso no haysuperficies características y en consecuencia las “zonas de Dependencia y de Influencia”coinciden y son iguales a todo el dominio del problema.Figura 4.- Zonas de influencia y de dependencia de un problema elíptico. La naturaleza del problema a resolver es un aspecto crucial en la selección de latécnica de discretización más apropiada. Hay que recurrir a la teoría de las ecuacionesdiferenciales en derivadas parciales. 27. Técnicas Numéricas en Mecánica de Fluidos 273.2.- Consideraciones físicas.Considérese un problema estacionario de transferencia de calor, en el que se requierehallar la distribución de temperatura en una barra de sección transversal rectangular constantey longitud infinita. El problema es bidimennsional y la ecuación a resolver es:∂2 T ∂2 T +=0 (20)∂ x2 ∂ y2 Para resolver el problema la ecuación debe ser completada con un conjunto decondiciones de contorno. Supóngase que se conoce una solución y que se introduce unaperturbación en un punto P (p.e. un foco de calor). Dicha perturbación va a ser percibida entodos los demás puntos de la barra. Por tanto, las temperaturas en todos los puntos deldominio están relacionadas entre sí y con las condiciones de contorno. Habrá que especificarla temperatura (o el flujo de calor) en todos los puntos del contorno para determinar lasolución. Se trata de un “problema de contorno”.Ejemplo: se quieren estudiar movimientos de pequeña amplitud en un conducto recto. Lasecuaciones gobernantes son:∂ρ∂ρ∂u +u+ρ=0∂t∂x∂x (21)∂u∂u1 ∂P +u=−∂t∂xρ∂xy la ecuación de la energía se reduce a la igualdad s = cte (entropía constante). Esto permiterelacionar las variaciones de presión y densidad:∂P∂ρ = a2(22)∂x∂xsiendo “a” la velocidad del sonido. Si la amplitud es pequeña: ∂ρ∂ρ∂u∂u u 0, mientras que adopta la forma: Fi +1φi +1 − Fi −1φi = Di +1 ( φi +1 − φi ) + Di −1 ( φi − φi −1 )(115)cuando uW < 0 y uE < 0. El mayor problema de esta discretización es el de su precisión cuando existanimportantes componentes de la velocidad en una dirección que no coincida con las líneas delmallado. El esquema aguas arriba produce entonces un “masajeado” artificial de las variacionesde las variables. El error resultante presenta la apariencia de una difusión y, por tanto,normalmente se le conoce como “falsa difusión”.Discretización híbrida.Se trata de una discretización que intenta solventar el problema de la discretización aguasF Farriba utilizando ésta cuando el cociente es mayor de 2 y una discretización central si< 2.D D 60. 60 Técnicas Numéricas en Mecánica de FluidosDiscretización según una ley polinómica.Desarrollada por Patankar en 1980. Constituye una metodología más precisa que laFhíbrida. No se crea falsa difusión si el cociente es superior a 10. En realidad sus propiedadesDson muy similares a las del algoritmo híbrido.Discretización QUICK. Se trata de una discretización de segundo orden que emplea una interpolacióncuadrática de tres puntos “aguas arriba” y ponderada para los valores de las funciones en lascaras de cada celda. En particular, la fórmula utilizada es:631 φi = φi −1 + φi − φi − 2 (116)888Tiene precisión de tercer orden para mallados uniformes y se trata de un algoritmoestable solo condicionalmente. Es decir, pueden aparecer problemas de pequeñas oscilacionesen la evolución de la solución.En cuanto a las técnicas para la evolución en el tiempo, es decir, técnicas explícitas oimplícitas (a veces aparece también el término pseudo-implícitas para describir una situaciónintermedia), se puede hablar también distintos métodos generales. Para una ecuacióndiferencial genérica donde la única variable sea la evolución temporal (variable “t”), es decir,del tipo: dU= S U + Q = H ( U, t )(117) dtdonde S, Q y, por tanto, H son funciones conocidas y U es la variable dependiente. Lasdiscretizaciones o métodos más utilizados son los siguientes: diferenciación hacia atrás, Eulerhacia atrás, trapezoidal de un paso y Euler explícito.De una forma genérica, las discretizaciones temporales se pueden expresar en función detres coeficientes base ξ, θ y ϕ para un instante “n” conocido, según la expresión:(1 + ξ ) U n + 2 − (1 + 2ξ ) U n +1 + ξU n = ∆ t [θ H n + 2 + (1 − ϕ − θ)H n +1 ] (118) Si se desea tener precisión de segundo orden en este tipo de discretizaciones, se debecumplir la relación específica:1ϕ = ξ−θ+(119)2A continuación se plantean las discretizaciones utilizadas según los distintos métodosmencionados. 61. Técnicas Numéricas en Mecánica de Fluidos 61Método hacia atrás (“backward”).Utiliza los valores ξ = 1/2, θ = 1 y ϕ = 0. Es decir, se cumple la ecuación 119. Suesquema general es, por lo tanto:3 n +22U 1 − 2 U n +1 + U n = ∆ t H n + 2 2[](120)Método Euler hacia atrás (“backward Euler”). En este caso, ξ = 0, θ = 1 y ϕ = 0. Es decir, no se cumple la ecuación 119. Su esquemageneral resulta con precisión de primer orden y es: (U n+2 − U n+1 ) = ∆ t [H n +2 ] (121)Método trapezoidal de un solo paso (también llamado de Crank-Nicolson). En este caso, ξ = 0, θ = 1/2 y ϕ = 0. Es decir, se cumple la ecuación 119 y es un esquemacon precisión de segundo orden. Su esquema general es: (U n+1 − U n ) = ∆2t [H n +1 + H n ] (122)Método Euler explícito. En este caso, ξ = 0, θ = 0 y ϕ = 0. Es decir, no se cumple la ecuación 119. Su esquemageneral resulta con precisión de primer orden y es: (U n+1 − U n ) = ∆ t [H n ](123)Ejemplo. Si se parte de la ecuación parabólica de difusión en dos dimensiones, es decir: ∂u ∂2 u ∂2 u = α 2 + ∂x (124) ∂t∂ y2  Si se discretiza en un mallado regular con ∆x = ∆y y utilizando un método centrado desegundo orden, para el término espacial, se llega a: ∂ u i, jα ∂t=∆x 2 ( u i+1, j + u i, j+1 + u i−1, j + u i, j−1 − 4 u i, j ) (125) 62. 62Técnicas Numéricas en Mecánica de FluidosAhora si se considera el método Euler hacia atrás para la discretización de la derivadatemporal, se obtiene:α ∆ t n +1 u i,+j1 = u i, j + n n ∆x 2( u i+1, j + u i,+j+11 + u in−+1,1 j + u i,+j−11 − 4 u i,+j1 ) n n n (126)Es decir, aparece un sistema penta-diagonal de ecuaciones lineales para cada nodo. Sifuera en coordenadas generalizadas (no ortogonales) sería un sistema con más de 5 elementosen la diagonal.o0o En realidad, si se aplica cualquiera de los métodos de evolución temporal y sedistingue entre las direcciones “x” e “y” (y si existe “z”) aparece un sistema tri-diagonal encada una de las direcciones, que se debe resolver por alguno de los métodos siguientes: • Métodos predictor-corrector. • Linealización local. • Métodos implícitos (ADI). • Métodos Runge-Kutta.4.6.- Resolución numérica de problemas sencillos. En este apartado se describen varias soluciones unidimensionales y bidimensionales decálculos realizados en problemas tipo, tales como la transmisión de calor en barras o el flujoen conductos.4.6.1.- Transferencia de calor en una barra alargada.Sea un problema de transferencia de calor, en el que se requiere hallar la distribuciónde temperatura en una barra de sección transversal rectangular y longitud mucho mayor quecualquier dimensión transversal, en la que se conocen las temperaturas en los dos extremos.El problema es unidimensional y la ecuación a resolver mediante el método de las diferenciasfinitas es: ∂2 T=0 (127) ∂ x2Para resolver el problema, la ecuación debe ser complementada con un conjunto decondiciones de contorno. Supóngase que se conoce una solución y que se introduce unaperturbación en un punto determinado (p.e. un foco de calor). Parece evidente que esaperturbación va a ser percibida en todos los demás puntos de la barra. Por tanto, lastemperaturas en todos los puntos del dominio están relacionados entre sí y con las condicionesde contorno. Habrá que especificar la temperatura (o el flujo de calor) en todos los puntos delcontorno para determinar la solución. Se trata de un problema de contorno. El ejemplo quese resolverá se muestra en la figura 26. 63. Técnicas Numéricas en Mecánica de Fluidos63 Figura 26.- Geometría de la barra unidimensional (21L >> L).Se lleva a cabo una discretización espacial considerando la barra como sólidounidimensional (tal y como se muestra en la figura 27). Figura 27.- Discretización para el problema unidireccional de transmisión de calor. Se toma como ejemplo práctico un caso con ∆x = 5 m y barra de aluminio(conductividad térmica k = 202.4 W/(mK)). Para el nodo i-ésimo se tendrá, utilizando unadiferencia centrada para discretizar la ecuación:Ti+1 − 2Ti + Ti-1=0 (128)2 ∆xEcuación que puede ser expresada en forma matricial: a1 a 20 0  T1   b1  T  b −1 2 1 0 2   2  0 −12 1 0  T3   b 3    0 ... 0    · =  (129)0 0   0 −1 2 −10        0 −1 2−1    b T20a 21   T   b 200 a 20  21   21  64. 64 Técnicas Numéricas en Mecánica de Fluidosdonde b2 = b3 = b4 = ... = b20 = 0 y los valores de a1, a2, a20, a21, b1 y b21 se deducen de lascondiciones de contorno (temperaturas en los extremos de la barra).Para los dos extremos se puede tomar una aproximación lineal o de segundo orden nocentrada (similar a la de la ecuación 128). Si se toma la aproximación lineal, se tiene (verfigura 28): T −T  T −T  k 2 1  = k 1 A (130) ∆x   ∆x/2  Figura 28.- Detalle de la discretización en la zona inicial de la barra.de donde se deduce: 2 TA − 3 T1 + T2 = 0 (131) De forma paralela, se tendrá:2 TB − 3 T21 + T20 = 0(132) Por lo que el sistema de ecuaciones completo será: 3 −100   T1   2 TA   −1 2 1 0  T2   0  0 −12 1 0  T3   0    0 ... 0    ·  =  (133)00     0 −1 2 −1 0       0 0 − 1 2 − 1   T0   20   T  0 − 1 3  21   2 TB  Conocidas las temperaturas en los extremos de la barra, se puede resolver el sistemade ecuaciones lineales descrito. En este caso, se obtiene una distribución lineal de temperatura 65. Técnicas Numéricas en Mecánica de Fluidos65entre A y B (independiente del material). Sean, por ejemplo, TA = 300 K y TB = 400 K, seobtiene:  T1   302.38      T2   307.14   T3   311.90         =(134)           T   390.48   20     T21   395.24  Esta distribución se representa gráficamente en la figura 29. T [K] 500 450 400 350 x [m] 300 0 204060 80 100 120Figura 29.- Distribución de temperatura para la conducción unidimensional.Si ahora, en vez de conocer TA y la temperatura en los dos extremos, se conoce unatemperatura y un flujo de calor, qB por ejemplo, la ecuación a resolver es la siguiente: ∂2 T k= q B . Dicha ecuación discretizada en el extremo B da lugar a: ∂ x2 T − T21 k  B  = qB (135) ∆x/2 es decir:qB TB = T21 −∆x(136) k 66. 66Técnicas Numéricas en Mecánica de Fluidosque introduciéndolo en la expresión de la ecuación para el último punto, da lugar al término: q − T20 + 3 T21 = 2 TB = 2  T21 − B ∆x  = 2 T21 − 0.00494 q B(137) k Si qB = 250 W/m2, se tendrá: − T20 + T21 = −1.23(138) Ahora se puede plantear de nuevo el sistema completo  T1   600  3 −1 0 0     T2   0  −1 2 1 0 T   0  0 −1 2 1 0  3    0 ...0·  =  (139)00    0 −1 2 −1 0    0 0 −12 −1       T20   0  0−1 1       T21   −1.23 que resuelto proporciona la solución mostrada en la figura 30. Se obtiene de nuevo unasolución lineal, pero con un mayor gradiente de temperaturas (TB = 430 K). T [K] 500 450 400 350 x [m] 300 0 20 4060 80100 120Figura 30.- Distribución de temperatura para barra unidimensionalcon flujo de calor en uno de sus extremos.Si existieran fuentes puntuales de calor, se pueden repetir los dos casos anteriores(temperatura en los dos extremo o temperatura en un extremo y condición de flujo de calor enel otro). La ecuación diferencial, para este nuevo problema, sería: 67. Técnicas Numéricas en Mecánica de Fluidos 67∂  ∂Tk +S = 0(140) ∂x  ∂x expresión en la que los coeficientes “S” serían los posibles valores de fuentes puntuales(términos independientes de las matrices).En concreto, se introduce una fuente puntual de calor en el nodo 11. Las ecuacionespara ambos casos no cambiarían y únicamente habría que introducir un valor distinto de 0 enel término independiente de la ecuación correspondiente al tramo número 11. Por ejemplo,para el caso de conocer las temperaturas en los dos extremos, se tendrá: 3 −1 0 0  T1   2TA   T   0 −1 21 0 2   0 −1  T3   0  2 1 0    0... 0   ·  =  b11  (141) 0 0   0 −1 2 −1 0         0 0 − 1 2 − 1       T20   0 0 − 1 3   T21   2TB  Si se particulariza para una fuente calorífica de valor S = 150 W/m3, el términoindependiente de la ecuación será: ∆x 5 mb11 = ( S ∆x )= ( 750 W/m 2 )= 18.53 K(142)k 202.4 W/(mK)con lo que el sistema a resolver es:  T1   2TA 3 -1 00     T2   0  -1 2 1 0  T   0 0 -1 2 1 0 3    0...0  ·   = 18.53  (143)0 0     0 -1 2 -1 0     0 0 -1 2 -1      T20   0  0 -1 3        T21   2TB Resolviendo para TA = 300 K y TB = 400 K, se obtiene la gráfica de evolución de latemperatura de la figura 31. 68. 68Técnicas Numéricas en Mecánica de Fluidos 4.60e+02 T [K] 4.40e+02 4.20e+02 4.00e+02 3.80e+02 3.60e+02 3.40e+02 3.20e+02x [m] 3.00e+02 Y0 20 40 60 80100120Figura 31.- Evolución de la temperatura en una barra con fuente de calor puntual. Esta distribución se puede mostrar también en forma de campo de temperatura a lolargo de la barra (figura 32), en la que se representa por medio de escala graduada latemperatura en la barra (escala de 300 K a 450 K). 4.50e+02 4.35e+02 4.20e+02 4.05e+02 3.90e+02 3.75e+02 3.60e+02 3.45e+02 3.30e+02 3.15e+02 YZ X 3.00e+02Figura 32.- Campo de temperaturas en una barra unidimensional con fuente puntual de calor.De forma totalmente paralela, el sistema para el caso de tener un flujo calorífico de150 W/m2 en el extremo B, resulta ser: 69. Técnicas Numéricas en Mecánica de Fluidos 69  3 −1 0 0  T1   2 TA       −1 21 0  T2   0   0 −1 2 1 0  T3   0        0... 0    0 0·  = 18.53 (144)      0 −1 2 −1 0           0 0 − 1 2 − 1     T0     20   0 − 1 1   T21   − 1.23    que para TA = 300 K, da lugar a la distribución de temperaturas mostrada en la figura 33,ahora con la escala de 300 K a 650 K. 6.50e+02 6.15e+02 5.80e+02 5.45e+02 5.10e+02 4.75e+02 4.40e+02 4.05e+02 3.70e+02 3.35e+02 YZ X 3.00e+02Figura 33.- Distribución de temperatura con flujo de calor en un extremo.Si existieran fuentes de calor no puntuales, no ocurriría como el caso anterior en el quesólo se altera la linealidad en un solo punto. Al introducir fuentes en más de un punto, ladistribución de temperaturas deja de ser lineal.Para los dos casos estudiados, se tendrán ecuaciones distintas. Si, por ejemplo, seintroducen fuentes puntuales en las 10 primeras celdas, se cambiarían esos diez primeroscoeficientes del término independiente (b1 a b10). De éstos, es distinto el correspondiente alprimer término, que se calcula según la figura 28, siguiendo la expresión:T −T  T −T k  A 1  − k  1 2  + q1 ∆x = 0(145) ∆x/2  ∆x es decir: 70. 70Técnicas Numéricas en Mecánica de Fluidoskk (2 TA − 2 T1 ) − (T1 − T2 ) + q1 ∆x = 0(146)∆x ∆xO, lo que es lo mismo: k k k k k3T1 = 2T1 +T2 + q1 ∆x = 2T1 +T + b1 (147)∆x∆x∆x∆x∆x 2Para el resto, se tiene, si se toma S1 = S2 =... = S10 = 150 W/m3:∆x 5 m b1 = b 2 = b3 = ... = b10 = ( S ∆x )= ( 750 W/m 2 )= 18.53 K (148) k 202.4 W/(mK) En el caso de que las temperaturas sean TA = 300 K y TB = 400 K, se obtiene elsiguiente sistema de ecuaciones, cuya solución se muestra gráficamente en la figura 34: 3 −1 00   T1   2 TA + 18.53    −1 21 0  T2   18.53  0 −1 2 1 0  T3   18.53    ... 0... 0  0 0·  =  18.53  (149)    00 −1 2 −1 0         ... 0 − 1 2 − 1    T0   20  0 − 1 3   T21     2 TB  9.00e+02T [K] 8.00e+02 7.00e+02 6.00e+02 5.00e+02 4.00e+02 x [m] 3.00e+02 Y020 40 60 80 100 120Figura 34.- Evolución de la temperatura en una barra con fuente de calor puntual. Esta distribución se puede mostrar también en forma de campo de temperatura a lolargo de la barra (figura 35), en la que se representa por medio de escala graduada latemperatura en la barra (escala de 300 K a 900 K). 71. Técnicas Numéricas en Mecánica de Fluidos71 9.00e+02 8.40e+02 7.80e+02 7.20e+02 6.60e+02 6.00e+02 5.40e+02 4.80e+02 4.20e+02 3.60e+02 YZ X 3.00e+02 Figura 35.- Campo de temperatura en una barra unidimensional con fuentes no puntuales de calor. En el caso que la temperatura sea TA = 300 K y el flujo de calor en B de 250 W/m2, elsistema a resolver es: 3 −10 0  T1   2 TA + 18.53   −1 2 1 0  T2   18.53  0 −12 1 0  T3   18.53   ...  0 ... 0    ·  =  18.53  (150)00    0 0 −1 2 −1 0   ... 0 − 1 2 − 1     T 0  20   T  0 − 1 1   21   − 1.23 La solución se muestra en forma de campo de temperatura a lo largo de la barra (figura36), en la que se representa por medio de escala graduada la temperatura en la barra (en estecaso, escala de 300 K a 1300 K).Al igual que en ejemplos anteriores, también se puede representar en gráficaunidimensional (figura 37). 72. 72 Técnicas Numéricas en Mecánica de Fluidos1.30e+031.20e+031.10e+031.00e+039.00e+028.00e+027.00e+026.00e+025.00e+024.00e+02 Y Z X3.00e+02Figura 36.- Campo de temperatura en una barra unidimensional con fuentes no puntuales de calor y flujo calorífico en un extremo. 1.30e+03T [K] 1.20e+03 1.10e+03 1.00e+03 9.00e+02 8.00e+02 7.00e+02 6.00e+02 5.00e+02 4.00e+02x [m] 3.00e+02Y0 20 4060 80100120Figura 37.- Evolución de la temperatura en una barra unidimensional con fuentes no puntuales de calor y flujo calorífico en un extremo.4.6.2.- Transferencia de calor en una placa plana bidimensional. En este caso, la ecuación a resolver es la de conducción de calor bidimensional en unaplaca plana, es decir:∂2 T ∂2 T +=0 (151)∂ x2 ∂ y2sobre la geometría de la figura 38 (se ha tomado una geometría rectangular con la base eldoble de la altura). 73. Técnicas Numéricas en Mecánica de Fluidos73 Figura 38.- Placa plana rectangular con relación de lados 0.5.Para el caso más elemental (sin fuentes de calor puntuales), haciendo un balance de losflujos que intervienen, se llega a: k i −1, j ∆ i −1, j( Ti, j − Ti−1, j ) − k ∆ ( Ti+1, j − Ti, j ) + i +1, j i +1, j x i, j − x i −1, j x i +1, j − x i, j    +  k i, j−1∆ i, j−1( Ti, j − Ti, j−1 ) − k ∆ ( Ti, j+1 − Ti, j ) =0 (152) i, j+1 i, j+1 yi, j − yi, j−1yi, j+1 − yi, j   donde los puntos (i-1, j), (i, j+1), (i, j-1) y (i+1, j) son los correspondientes a la nomenclaturade la figura 39.Figura 39.- Celda unitaria para la discretización bidimensional de un dominio. 74. 74 Técnicas Numéricas en Mecánica de Fluidos Si k es constante y se toma un mallado uniforme, es decir: ∆i-1, j = ∆i+1, j = ∆i, j-1 = ∆i, j+1= ∆x = ∆y = ∆, se obtiene:k  4 Ti, j − Ti +1, j − Ti −1, j − Ti, j+1 − Ti, j−1  = 0 (153)de donde: 4 Ti, j = Ti +1, j + Ti −1, j + Ti, j+1 + Ti, j−1  (154) Fórmula que generalizada para cualquier punto “i, j”, proporciona la expresión:Ti j = 1 [T + Ti+1 j + Ti j−1 + Ti j+1 ] 4 i −1 j (155) Como condición de contorno, se impondrá la temperatura en los laterales de la placaplana. En particular se supondrá que la temperatura es uniforme en todo el contorno exceptoen una zona central situada en la base de la placa, tal y como se muestra en la figura 40. Figura 40.- Condición de contorno para el problema de la conducción en una placa plana.En concreto, se resuelve el problema para obtener la distribución de temperaturas enuna placa de (6×3 m2, es decir, L = 3 m) y considerando una discretización de 24×12 celdas.Se impone una zona central formada por 8 celdas (2 m) en la base con una mayor temperaturaque el resto. Las condiciones de contorno elegidas son: 75. Técnicas Numéricas en Mecánica de Fluidos75 T0 j = 300 K ∀ j Ti12 = 300 K ∀ i T24 j = 300 K ∀ j (156) Ti 0 = 300 K para i = 0...7 e i = 17...24 Ti 0 = 800 K para i = 8...16 Se ha resuelto el problema por dos procedimientos distintos: resolución directa delsistema de ecuaciones lineales y utilizando un código numérico comercial. Las dos solucionesobtenidas para este problema se muestra en las figuras 41 y 42 para distintas posiciones “i”. En la figura 41 se muestra la solución directa del sistema lineal de ecuaciones. 750 T [K] y=2.75 700 y=2.5 m 650 Y=2.25 m 600 Y=2 m 550 y=1.75 m 500 y=1.5 m 450 y=1.25 m y=1 m 400 y=0.75 m 350 y=0.5 m 300 y=0.25 m 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23Posiciones “i” Figura 41.- Conducción de calor en una placa plana. Distribución de temperaturas. En la figura 42 se muestra la solución obtenida mediante un código numéricocomercial. Además, con este código se puede obtener una distribución bidimensional delcampo de temperaturas en la placa. Para las condiciones de contorno de la ecuación 156, seobtiene el resultado mostrado en la figura 43. 76. 76 Técnicas Numéricas en Mecánica de FluidosT [K]x [m]Figura 42.- Conducción de calor en una placa plana. Distribución de temperaturas.Figura 43.- Campo de temperaturas en el problema de conducción en placa plana.Para la solución obtenida mediante el código numérico se obtiene (comparar figuras41 y 42) una distribución con un mayor gradiente de temperaturas en la zona de calentamientoen la base, debido a que en este caso en la discretización espacial se han utilizado más celdas. Una segunda condición de contorno probada es la de imponer un flujo de calor en lasparedes verticales y en la pared superior. En concreto se elige la condición de flujo nulo, queequivale a un aislamiento externo adiabático.Para esta nueva condición, una vez discretizado y considerando la ecuación queresultaría para la pared derecha (lado “i-1” en la figura 39), se obtiene:2  Ti −1, j Ti, j−1 Ti, j =  + (157)3  22  77. Técnicas Numéricas en Mecánica de Fluidos77Impuesta esta nueva relación en el sistema a resolver, se llega a la distribución detemperaturas mostrada en la figura 44. 800 T [K]y=2.75 750y=2.5 m 700Y=2.25 m 650Y=2 m 600y=1.75 m 550y=1.5 m 500y=1.25 m 450y=1 m 400 350y=0.75 m 300y=0.5 m 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23y=0.25 m Posiciones “i” Figura 44.- Conducción de calor en una placa plana conflujo neto nulo en tres paredes. Distribución de temperaturas.Los resultados que se obtienen con el código numérico se muestran (ver figura 45) enforma similares distribuciones a las mostradas en la figura 44 y en forma de campo detemperaturas en la figura 46. T [K]x [m] Figura 45.- Conducción de calor en una placa plana con flujo neto nulo en las paredes. Distribución de temperaturas obtenida utilizando el código numérico. 78. 78Técnicas Numéricas en Mecánica de Fluidos Figura 46.- Campo de temperaturas en una placa planacon condiciones de flujo nulo en tres caras.La solución a partir del sistema de ecuaciones se ha implementado en una hoja decálculo en forma de algoritmo que resuelve el caso de la transmisión de calor en una placabidimensional tanto imponiendo condiciones de temperatura, como imponiendo condicionesde flujo de calor. En la figura 47 se muestra la salida que proporciona dicho algoritmo en elúltimo caso. 12 300 300 300 300 300 300 300 300 300 300 300 300 300 300 300 300 300 300 300 300 300 300 300 300 300 11 300 302 304 306 309 311 313 315 317 319 321 322 322 322 321 319 317 315 313 311 309 306 304 302 300 10 300 304 308 312 317 322 326 331 335 339 342 344 344 344 342 339 335 331 326 322 317 312 308 304 300 9300 306 312 318 325 332 340 347 354 360 364 367 368 367 364 360 354 347 340 332 325 318 312 306 300 8300 308 316 324 333 343 353 363 373 381 388 392 393 392 388 381 373 363 353 343 333 324 316 308 300 7300 309 319 329 341 353 366 380 394 405 414 420 422 420 414 405 394 380 366 353 341 329 319 309 300 6300 310 321 333 347 362 379 398 416 432 444 451 454 451 444 432 416 398 379 362 347 333 321 310 300 5300 311 322 336 351 370 391 415 440 462 478 488 491 488 478 462 440 415 391 370 351 336 322 311 300 4300 311 322 336 352 374 401 433 467 497 519 531 535 531 519 497 467 433 401 374 352 336 322 311 300 3300 309 320 333 349 372 405 448 499 541 568 583 588 583 568 541 499 448 405 372 349 333 320 309 300 2300 307 315 325 340 362 397 455 542 598 629 645 650 645 629 598 542 455 397 362 340 325 315 307 300 1300 304 308 314 323 338 367 434 614 680 707 718 722 718 707 680 614 434 367 338 323 314 308 304 300 0300 300 300 300 300 300 300 300 800 800 800 800 800 800 800 800 800 300 300 300 300 300 300 300 30001 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24Figura 47.- Solución del sistema de ecuaciones en una discretización 23×11. 79. Técnicas Numéricas en Mecánica de Fluidos 794.6.3.- Flujo entre placas planas.Se trata aquí de resolver el problema del flujo de un fluido entre dos placas planas,siendo la separación entre ambas mucho menor que la longitud de los lados de dichas placas(figura 48). yh xLFigura 48.- Geometría para el problema del flujo entre placas planas paralelas. Se parte de las ecuaciones de Navier-Stokes (ecuación 2) y se introducen las siguienteshipótesis simplificativas: i) Flujo incompresible: ∇·u = 0 . ii) Flujo estacionario, variaciones temporales nulas. iii) Movimiento bidimensional (separación entre placas mucho menor que longitud delas mismas, h 0 o bien en la figura 51) ocontrario (p* < 0 o bien en la figura 52) que el movimiento de la placa.yUUh + = xL Figura 51.- Perfil parabólico más perfil lineal para el flujo entre dos placas con p* > 0.yUUh + = x L Figura 52.- Perfil parabólico más perfil lineal para el flujo entre dos placas con p* < 0. 82. 82 Técnicas Numéricas en Mecánica de Fluidos Se puede obtener el caudal circulante integrando la velocidad entre las dos placas: Q h U h p* h 3q= b = ∫ 0 u(y) dy =-2 12 µ (170)que se anula (q = 0) para la condición de p* negativo e igual a:6µ U p* = −(171)h2A continuación, se particulariza para la geometría mostrada en la figura 53 (h = 0.005m y L = 5 m) con ancho b constante (perpendicular al plano de dicha figura) y para un fluidocon las propiedades siguientes: ρ = 889.2 kg / m3 y µ = 1·10−6 Pa s .5 metros 5 mm Figura 53.- Solución numérica particular para el flujo entre placas paralelas.Además de haber obtenido las expresiones teóricas exactas, se puede resolver el flujoutilizando un método numérico (código comercial Fluent). Para el caso de flujo de Couette yvelocidad uniforme de 5 m/s y se obtiene el resultado mostrado en la figura 54.Figura 54.- Flujo con perfil lineal entre dos placas planas. Para el caso de flujo de Pouseuille con ∆p = 8.5 MPa (p* = ∆p/L = 1.7 MPa/m), seobtiene el perfil parabólico de la figura 55. Se ha escogido este valor porque según laecuación 168 produce una velocidad máxima de 5 m/s. 83. Técnicas Numéricas en Mecánica de Fluidos83Figura 55.- Flujo parabólico entre dos placas planas paralelas. Si se resuelve la situación mostrada en la figura 51, es decir, gradiente de presiónsuperpuesto al perfil lineal y produciendo flujo el mismo sentido, se obtiene la soluciónmostrada en la figura 56. Resolviendo la ecuación 169 se obtiene una velocidad máxima de7.82 m/s, que es prácticamente la misma que la que proporciona el método numérico.Figura 56.- Perfil parabólico superpuesto al lineal con efectos en el mismo sentido.En la figura 57 se muestra la distribución de velocidad para las distintas alturas y entrelas placas (0 < y < 0.005 m), observándose claramente la posición y valor del máximo. 84. 84Técnicas Numéricas en Mecánica de Fluidosu (y) [m/s]y [m]Figura 57.- Distribución de velocidad entre las dos placas planas con superposición de efectos.Si, por el contrario, el efecto del gradiente de presión es opuesto al del movimiento dela placa superior (ecuación 169 con p* < 0), se obtiene el resultado de la figura 58. Se utilizanlos mismos valores numéricos: velocidad uniforme de 5 m/s y ∆p = 8.5 MPa.Figura 58.- Perfil parabólico superpuesto al lineal con efectos en sentido opuesto. En la figura 59 se muestra la evolución de la velocidad para una posición “x” dada enfunción de la altura. 85. Técnicas Numéricas en Mecánica de Fluidos85 u (y) [m/s]y [m] Figura 59.- Distribución de velocidad entre las dosplacas planas con superposición de efectos opuestos. Se puede resolver este caso para la condición de caudal nulo (Q = 0), en cuyo caso elgradiente de presión reducida ha de ser:6µ U p* = −= 1.272 MPa (172) h2 Imponiendo un ∆p = 6.36 MPa, se obtiene la distribución de velocidad y el perfil enuna posición “x” mostrados en las figuras 60 y 61. Figura 60.- Distribución de velocidad para el caso de Q = 0. 86. 86Técnicas Numéricas en Mecánica de Fluidosu (y) [m/s] y [m] Figura 61.- Perfil de velocidad en una posición “x” en función de la altura “y”.Numéricamente, a partir de la solución obtenida, se comprueba que Q = 0 para el casoasí definido.4.6.4.- Flujo en un conducto cilíndrico. Se resuelve ahora el flujo laminar en un conducto cilíndrico. A partir de las hipótesisya señaladas para el flujo entre placas planas (flujo estacionario, incompresible) y añadiendola condición de axi-simetría impuesta por la geometría cilíndrica, se puede obtener, para latensión cortante:r τ = − p* (173)2donde, al igual que en el apartado anterior, p* es el gradiente de presión reducido, que nodepende del radio considerado y que se puede expresar en función de la nomenclatura de lafigura 62 como: ∆pp* =+ ρ g cos θ (174) ∆L Expresión que tiene en cuenta el efecto gravitatorio y la diferencia de presiones en elconducto entre las dos secciones “1” y “2” separadas una distancia ∆L. Se obtiene la siguiente distribución de velocidad: u= 4µ( p* 2r − R2 )(175) 87. Técnicas Numéricas en Mecánica de Fluidos 87 r p1 p2θ D ∆L Figura 62.- Geometría básica para el flujo laminar en conductos cilíndricos.que constituye un perfil parabólico y que, habitualmente, se expresa en función de la velocidadmáxima, según la expresión:  R 2p*R2 u = u MÁX 1 −    con u MÁX = − (176)  r 4µ El caudal circulante será entonces:πp* 4 Q=− R (177)8µ A continuación se va a resolver este problema mediante un código comercial. Se resuelveel caso particular de tener un fluido con unas propiedades dadas por:ρ = 889.2 kg/m3 y µ = 1·10-6 Pa s .y un conducto horizontal, tal y como se muestra en la figura 63.5 metrosØ 5 cmFigura 63.- Geometría para el problema del flujo en conductos cilíndricos. Para ello se discretiza el dominio usando un mallado no estructurado en la seccióntransversal y proyectándolo en la dirección del flujo, tal y como se muestra en la figura 64. 88. 88 Técnicas Numéricas en Mecánica de Fluidos Figura 64.- Geometría discretizada para resolución del caso de conducto horizontal. El movimiento del flujo será debido únicamente a la diferencia de presión entre laentrada y la salida. Se fija un caudal que asegure el flujo laminar, por ejemplo, 3.92 l/s. Elgradiente de presión reducida necesario será: 8µ Qp* == 1693 Pa / m(178) πR4 La velocidad máxima, que se alcanza en el eje del conducto, será:p*R 2u MÁX = − = 1 m/s(179) 4µLa resolución del flujo en el interior del conducto permite visualizar la velocidad endistintas secciones (figura 65). En dicha figura se muestran los vectores velocidad en seissecciones distintas del conducto. Figura 65.- Distribución de velocidad en un conducto circular con flujo laminar. 89. Técnicas Numéricas en Mecánica de Fluidos 89 Tal y como se había previsto con la ecuación teórica, se obtiene un perfil parabolico develocidad. En la figura 66 se muestra dicho perfil en la sección de salida, donde cada puntorepresenta el valor de la velocidad en cada celda de la sección transversal del conducto. u (r) [m/s] r [m]Figura 66.- Distribución de velocidad en la salida (perfil parabólico). También se puede visualizar la presión en distintas secciones, observándose sucarácter lineal con la distancia (figura 67).Figura 67.- Distribución lineal del salto de presión. ∆P= 8690 Pa. En caso de utilizar un conducto con inclinación, la diferencia de presión varía alintroducirse el término ρg cosθ en la ecuación 174, pero el caudal, que es el que impone unavelocidad máxima del flujo, y el perfil parabólico se mantienen.Con los datos de partida se obtiene un número de Reynolds de 84. En caso de tener unnúmero de Reynolds mayor de 4000, el flujo deja de ser laminar y se hace necesaria lautilización de un modelo de la turbulencia para la resolución del flujo. 90. 90Técnicas Numéricas en Mecánica de Fluidos4.6.5.- Flujo incompresible no estacionario: apertura de válvula en conducto. Los flujos incompresibles han sido mencionados como una aplicación particular de lasecuaciones de Navier-Stokes en apartados anteriores. Se profundiza un poco más en suestudio y se presenta un ejemplo de cálculo basado en el algoritmo del paso fraccional.A continuación se describe la resolución de uno de los problemas clásicos en laMecánica de Fluidos: la descarga de un depósito a través de un conducto con una válvula quetiene una apertura dependiente del tiempo. Inicialmente se describe la problemática general delos flujos incompresibles para pasar luego a la resolución del problema concreto.Introducción: Estudios sobre flujo incompresible mediante métodos numéricos.En el artículo de Gresho (1991), publicado en el Annual Review of Fluid Mechanics,se realiza una exposición muy detallada sobre las posibles formas de abordar numéricamenteun problema de flujo incompresible. Primero expone la formulación matemática genérica de este tipo de problemas parapasar luego a señalar las opciones en cuanto al planteamiento de los distintos términospresentes en las ecuaciones de Navier-Stokes. Distingue así las siguientes formulaciones: lade divergencia de las tensiones, la clásica, la mixta divergencia-rotacional, la rotacional, laconvectiva difusiva y la simétrica, entre otras. Añade a estas formulaciones básicas dos más,derivadas de las anteriores, con gran aplicación práctica, como son la formulación de laecuación de Poisson para la presión (PPE) y la de la ecuación de transporte para la vorticidad(VTE).Después discute las distintas posibilidades que conducen al establecimiento de unproblema “bien definido” en lo que se refiere a las condiciones de contorno y las condicionesiniciales. Particulariza todas las posibilidades para los tipos de formulación de las ecuacionesdiscutidas anteriormente. Finalmente, realiza algún estudio más a fondo sobre distintosproblemas tipo, particularizando para problemas en los que hay saltos bruscos en lascondiciones de contorno.Se trata de un artículo básico en el estudio numérico de los flujos incompresibles y, sibien, su orientación es hacia los flujos no estacionarios, la mayor parte de las conclusiones estotalmente aplicable a flujos estacionarios, considerando a estos como una evolución noestacionaria que se estabiliza en un valor constante. Kwak en sus múltiples artículos (Kwak et al., 1986; Kwak et al., 1992; Rogers et al.,1989 o Rosenfeld et al., 1991) y, sobre todo, en el monográfico sobre métodos numéricosusados en el cálculo de flujos incompresibles (Kwak, 1989) analiza un amplio rango deposibles aplicaciones para el estudio de flujos bi y tridimensionales. Ha desarrollado numerosos algoritmos basados en formulaciones que emplean lasvariables primitivas (presión y velocidad). Destacan la aplicación práctica de algoritmos depseudocompresibilidad (Kwak et al., 1986) y de paso fraccional (Rosenfeld et al., 1991). Ensu trabajo, ha presentado numerosas validaciones con problemas tipo (entre otros, flujo en unescalón, flujo en una cavidad cuadrada y flujo en la estela de un cilindro), así como multitud 91. Técnicas Numéricas en Mecánica de Fluidos91de soluciones numéricas a problemas prácticos tan dispares como el flujo en el motor de unalanzadera espacial o el flujo en un implante de un corazón. Constituye un referente ineludible en cualquier estudio relacionado con la aplicaciónde técnicas numéricas al estudio de flujos incompresibles.Posibilidades de resolución de las ecuaciones para flujos incompresibles no estacionarios.El desarrollo de las técnicas numéricas ha permitido un análisis cada vez más rigurosode los flujos, mejorando principalmente lo que se refiere a flujos compresibles. Los flujosincompresibles han seguido esta evolución con cierto retraso, aún cuando muchas de lasrealizaciones industriales son aplicación de flujo incompresible. Dos dificultades básicas hanfrenado el desarrollo de los métodos numéricos aplicados a flujos incompresibles respecto a losequivalentes para flujo compresible. Por un lado la excesiva y ya descrita “rigidez” de lasecuaciones y, por otro, el menor interés de estos flujos en el campo de la aeronáutica militar.Así, la mayor parte del esfuerzo investigador se ha dedicado a modelizar el estado estacionariode la evolución del fluido en los distintos equipos, mientras que muchos flujos reales sonesencialmente no estacionarios (desprendimiento de vórtices, turbulencia, flujo enturbomáquinas o en sistemas biológicos, entre otros).Tal y como señalan Peyret y Taylor (1983) o Hirsch (1988), existen dos grandes gruposde métodos según se plantee la resolución de las ecuaciones del flujo incompresible: las basadasen el empleo de variables primitivas y las que emplean la ecuación de la vorticidad. Estasúltimas surgieron debido a la inexistencia de condiciones de contorno físicas de contacto consólidos para la presión. Sin embargo, la búsqueda de una condición de contorno para lavorticidad, salvo en casos puntuales, resulta más complicado.A pesar del gran abanico de variantes que presentan ambas formas de solución (verQuarteroni, 1996), son básicamente las que usan variables primitivas las únicas usadas para elcaso de flujo no estacionario. Entre éstas destacan los algoritmos que operan con la ecuación dePoisson para resolver el campo de presiones, y el método de pseudocompresibilidad. Ambostratan de solventar la “rigidez” de las ecuaciones para flujo compresible a bajos números deMach. La diferencia básica entre los dos métodos reside en el carácter de las ecuaciones aresolver.En el caso de la compresibilidad artificial o pseudocompresibilidad se trata de adaptar alcálculo de flujos incompresibles las soluciones que se tienen para el caso de fluidoscompresibles. El método introduce un parámetro de compresibilidad artificial para lograr laadaptación de las ecuaciones (Kwak, 1992). El sistema resultante se resuelve utilizando unalgoritmo de “evolución en el tiempo” (time-marching). Su principal aplicación ha sido lamodelización de casos estacionarios, debido a que sólo mediante modificaciones en el saltotemporal, introduciendo tiempos artificiales (pseudotime) o “inclinación temporal” (time-inclining), se ha podido llegar a soluciones no estacionarias. La inclusión en el modelo de unparámetro artificial, y la crucial dependencia de algunas características del algoritmo numéricocomo son la convergencia o la estabilidad con dicho parámetro, hacen que el método se puedaconsiderar poco adaptable a problemas no estacionarios. 92. 92 Técnicas Numéricas en Mecánica de FluidosRespecto a las soluciones que emplea la ecuación de Poisson para la presión, se puedendistinguir dos aproximaciones diferentes, de características similares, pero de muy distintaimplementación:i) Método de presión corregida. Introducido por Patankar (1980), consiste en unamanipulación de la presión que da lugar a un proceso de solución iterativo para los campos depresiones y velocidad. La aplicación de este razonamiento, ha dado lugar a un conjunto desoluciones para el caso de flujos estacionarios usado en muchos programas comerciales. Para elcaso no estacionario su mayor desventaja es la necesidad de condiciones de contorno para elcampo de presiones ficticio y la de requerir una iteración adicional para cada paso temporal. Aeste tipo de métodos se les denomina SIMPLE, debido a que son estas las iniciales delalgoritmo más importante desarrollado por Patankar (1980). ii) Método de paso fraccional. Desarrollado a partir de los algoritmos MAC (Mark andCell), su aplicación para el estudio de flujos incompresibles y no estacionarios fue establecidapor Peyret y Taylor (1983). Se parte de la formulación de ecuación de Poisson para la presión yse trata de discretizar las ecuaciones de gobierno para luego diferenciar la ecuación de cantidadde movimiento aplicando la condición de divergencia nula para la velocidad. El método de pasofraccional es una de las posibilidades más aconsejadas en la bibliografía consultada (Rosenfeldet al., 1991 o Gresho, 1991) para el tratamiento de flujos no estacionarios. Además, a través deestrategias computacionales como el multimallado, puede conducir a soluciones con un tiempode cálculo muy reducido, lo cual le daría una gran ventaja respecto al método de lapseudocompresibilidad.Ejemplo de algoritmo de resolución. Se considera primeramente la simplificación de fluido ideal (µ = 0) y se plantea unaresolución de las ecuaciones de gobierno para flujo incompresible siguiendo un modelo de pasofraccional. Las ecuaciones de Euler de continuidad y de cantidad de movimiento para flujoincompresible son:∇·u = 0(180)Du 1 =-∇p+g(181)Dt FrEcuación en la que se han usado variables adimensionales y p representa la presióncinemática, es decir, presión dividida por densidad. El desarrollo del método se va a realizarconsiderando que las fuerzas gravitatorias son despreciables. Discretizando la ecuación decantidad de movimiento entre dos pasos temporales consecutivos, “n” y “n+1”, y usandoalgoritmo explícito con primer orden de aproximación temporal se obtiene: n +1 -u n u+ ∇ p n +1 + ( u n · ∇ ) u n = 0 (182)∆t Se introduce ahora una velocidad auxiliar intermedia, denominada ficticia ( u f ), deforma que se cumple la relación siguiente: 93. Técnicas Numéricas en Mecánica de Fluidos93 fnu - u + ( n · ∇) n = 0(183) uu ∆t Según la ecuación de cantidad de movimiento 182, esta última relación obtenidaimplica:n +1u-uf + ∇ p n +1 = 0 (184)∆t La condición que se busca para el paso temporal “n+1” es la adivergencia del campo develocidades. Expresada matemáticamente, dicha condición se escribe: ∇· u n +1 =0 (185)La forma de operar del método de paso fraccional se puede resumir en los siguientescuatro pasos básicos: i) Resolución de la ecuación 183, obteniéndose de forma explícita el campo de velocidades intermedio ( u f ), por ser la única incógnita en el paso temporal n. ii) Se toma la divergencia de la ecuación 184 e imponiendo la condición dada por 185, se obtiene la ecuación de Poisson para la presión: 1∆ p n +1 =∇·u f (186) ∆t iii) Una vez resuelta esta última, se puede obtener la solución para el campo de presiones. iv) Resolviendo 184, con pn+1 ya conocido, se llegaría al campo de velocidades adivergente buscado para paso temporal n+1.Algoritmo unidimensional. Problema de la descarga no estacionaria de un depósito.La descarga de un depósito constituye un problema clásico en Mecánica de Fluidos y hasido tratada ampliamente en la bibliografía técnica (Liñán, 1991 o Pentaris et al. 1996). A pesarde la aparente sencillez, se trata de un problema con escalas, tanto espaciales como temporalesmuy diferentes, lo que posibilita distintas aproximaciones y a la vez imposibilita una soluciónexacta. Aquí se resolverá para una condición de contorno no estacionaria consistente en uncierre periódico de la válvula de salida.Usando una discretización según diferencias finitas para el problema unidimensional, sepuede llegar a una discretización para la ecuación 186. Se trata de un sistema de ecuacionestridiagonal con incógnitas las presiones en tres puntos (a indica los valores de las áreas del canalen los distintos puntos). Dicho sistema tridiagonal (ecuación 137) se resuelve numéricamente,obteniéndose el campo de presiones para cada paso temporal considerado (paso “n + 1”). 94. 94 Técnicas Numéricas en Mecánica de Fluidos∆xn+1n+1n+1 A i-1/2 pi-1 - (A i-1/2 + A i+1/2) pi + A i+1/2 pi+1 =  A i+1/2 u fi+1/2 - A i-1/2 u fi-1/2   (187)∆t Con el campo de presiones calculado, se obtiene la velocidad en cada punto y enconcreto para el punto de salida del depósito. El problema así planteado posee una solución analítica que se expresa en función de lafrecuencia de cierre de la válvula ω y de la velocidad adimensional en la sección de salida deldepósito hacia el conducto (u), que se obtiene resolviendo la ecuación diferencial: ln ( AR ) d u 11 + 2u = [ 1 − cos (ω t)](188) AR − 1 d t 2 AR24 Donde “AR” es la relación de anchos entre la secciones de entrada y de salida del tubode descarga. Se observa como detalle importante que la ecuación diferencial del fenómeno no seve afectada por la longitud total del canal. Este hecho es coherente con haber considerado lasolución de las ecuaciones de Euler. En la figura 68 se puede comprobar la comparación entre la solución analítica y lasolución usando el modelo y discretización explicados. Ambas evoluciones se muestran pormedio de la variable u*, que sería equivalente a la variable u en la ecuación 187, en funcióndel tiempo adimensionalizado (t*), también equivalente a la variable t en 187. Tal y comopuede verse, la concordancia de ambos resultados es muy grande.1.00 u*0.750.50 Modelo Unidimensional0.25 Ecuación Diferencial0.0005 10 15 t* 20Figura 68.- Comparación de resultados entre la resolución teórica (ecuación 140) y el modelo unidimensional (González et al., 1997). 95. Técnicas Numéricas en Mecánica de Fluidos 95Algoritmo bidimensional. Según lo dicho al comienzo de este apartado, es posible resolver un flujoincompresible usando únicamente las ecuaciones de continuidad y cantidad de movimiento.Como modelo físico, en este apartado se plantean las ecuaciones de gobierno para flujobidimensional, incompresible no estacionario, de un fluido ideal, expresadas en un sistema dereferencia curvilíneo ortogonal (ξ,η). A través del cambio de variable se posibilita resolver elproblema utilizando una técnica de diferencias finitas. De forma adimensional y para el caso de fuerzas másicas nulas, se tiene:  ∂ ( )  ∂ ξ h2 u + ∂ η∂ ( h v) = 0 1   1  ∂ ∂p  ∂ u  ∂t +  h1 h2  ∂ ξ (h u)+22∂∂η ( h u v) + u v ∂∂ hη 11− v2 ∂ h2 ∂ξ+ h2 = 0 ∂ξ     ∂ v1  ∂ ∂p   +  ( h u v) +1 ∂ ( h v ) + u v ∂∂ hξ122− u2 ∂ h1+ h1 = 0 ∂ th1 h2  ∂ ξ ∂η ∂η∂η (189) Ecuaciones en las que u y v son las componentes de la velocidad ( u ) en función deltiempo (t), expresadas en la referencia del sistema ortogonal y los hi son los coeficientes métricoscorrespondientes al cambio de base del sistema de coordenadas, necesario para poder continuarusando la discretización en diferencias finitas. En este caso, se realiza una discretización espacialsegún diferencias finitas que asegure una aproximación de segundo orden al sistema deecuaciones diferenciales.Se construye un mallado bidimensional escalonado (staggered grid), que se representaen la figura 69. Los términos espaciales se discretizan mediante diferencias centradas, exceptoen los puntos del contorno. Figura 69.- Esquema del mallado escalonado utilizado. 96. 96 Técnicas Numéricas en Mecánica de Fluidos Según la nomenclatura de la figura 69, la discretización de la ecuación general para laprimera componente de la velocidad resulta ser:u if+1 / 2, j = u in+1 / 2, j − ∆th1i +1 / 2, j ∂ [ h2un (h 2 i +1 / 2, j  ∂ ξ i +1 / 2, j i+1/ 2 , j ) ]+ 2∂h1 (u n v in+1 / 2, j ∂ η i +1 / 2, j i+1 / 2 , j)  ∂ h1i +1 / 2, j   2 ∂ h2 i +1 / 2 , j  + u n+1 / 2 , j v in+1 / 2, j    − (v in+1 / 2, j )    i ∂η  ∂ξ       (190) Tras resolver esta ecuación y la equivalente para la segunda componente de la velocidadf(v ), se plantea la ecuación de Poisson para la presión, que será en este caso: ∂  h2 i,j ∂ pi , j  h + ∂  1 i , j  ∂ pi , j   = 1  ∂  h 2∂     u fi +1 / 2, j  + h1 i , j+1 / 2 u fi, j+1 / 2  ∂ ξ  h1 i , j  ∂ ξ  ∂ η  h ∂ η ∆ t  ∂ ξ  i+1 / 2 , j  ∂ η    2 i,j  (191)Tras discretizar las derivadas planteadas, se llega a un sistema penta-diagonal deecuaciones lineales. Dicho sistema se resuelve usando un método aproximado, a diferenciadel sistema tridiagonal que se obtuvo en el caso unidimensional, mostrado anteriormente. Seelige un método implícito según direcciones alternativas con aceleración de la convergencia(Accelerated Alternating Direction Implicit, conocido por AADI). La aceleración en laconvergencia se hace por medio de un parámetro λ de sobrerrelajación, que reducenotablemente el número de iteraciones en cada paso temporal. En la figura 70 se presenta ladisminución tanto en tiempo de cálculo como en número de iteraciones para uno de los casosprobados y la forma de obtener el parámetro λ óptimo, que resulta ser función de ladiscretización espacial, y más en concreto del parámetro β (β=∆ξ/∆η). λ1.41.3 λopt=1.2711.21.11.0Iteraciones Tiempo de CPU [s]0.90.8 0 50 100150200250300Figura 70.- Evolución del número de iteraciones y del tiempo de cálculoen función del parámetro de sobrerrelajación (λ). 97. Técnicas Numéricas en Mecánica de Fluidos97Una vez resuelto el campo de presiones en cada instante temporal, se pasa a actualizar elcampo de velocidades como segundo paso del método de proyección (se proyecta la velocidadhasta conseguir su adivergencia). La discretización para la obtención de la primera componentede la velocidad resulta:+ u in+11/ 2, j = u f i +1 / 2, j −∆t 2 ∆ξ h 1i +1 / 2, j (++ p in+11j − p in−11j,, ) (192)Con esta última operación concluye el cálculo para cada instante. La precisión delmétodo descrito es de primer orden para la discretización temporal y de segundo orden en ladiscretización espacial, es decir, los residuos son del orden de O(∆ξ2, ∆η2, ∆t).Se obtienen también los correspondientes campos de presión y velocidades en cadainstante para la geometría considerada (difusor recto). La solución para distintas relaciones deárea del difusor puede verse en la figura 71. Dada la independencia del modelo del ángulo deldifusor (modelo no viscoso), se han representado los campos de velocidades para unageometría del difusor con ángulo θ = 30o, aunque para las relaciones de área estudiadas ellímite práctico en el ángulo del difusor es mucho menor (2θ ≈ 15o). 1.6 u* AR = 2.0(∆t=2·10E-4) 1.4 1.2AR = 1.5(∆t=2·10E-3) 1.0 0.8 AR = 1.0 (∆t=1·10E-2) 0.6 0.4 Ec. Dif. 0.2 Modelo 2D t* 0.0 0 5 10 152025Figura 71.- Comparación de resultados entre la resolución teórica y el modelo bidimensional (González et al., 1998). 98. 98 Técnicas Numéricas en Mecánica de Fluidos 99. Técnicas Numéricas en Mecánica de Fluidos 995.- Discretización del dominio: generaciónde mallados.Las técnicas numéricas para la generación de un mallado sobre el que resolver lasecuaciones de gobierno se definen (Thomson et al., 1982; Thomson et al., 1983; Thomson,1984 o Niederdrenk, 1987) como procedimientos para la distribución ordenada deobservadores o estaciones de toma de datos en un dominio físico, de forma que exista laposibilidad de comunicación eficiente entre los distintos observadores y que los distintosfenómenos físicos que ocurran en dicho dominio físico continuo estén representados con lasuficiente exactitud por medio de dichos observadores (discretización espacial).Los mallados se pueden clasificar atendiendo a dos criterios: la forma de definir lasfronteras del dominio y la conectividad entre los distintos observadores o puntos del mallado. En cuanto a la forma de definir las fronteras, están los mallados conformes (nointerpolan en la definición de la frontera) y los no conformes (que no definen la frontera realsino una aproximación numérica de la misma).Atendiendo al criterio de conectividad o estructura de los datos del mallado, existendos métodos de generación de mallados que darían lugar a los dos tipos de mallados básicos:los estructurados y los no estructurados. En los mallados estructurados, los observadores secolocan siguiendo una red de familias de líneas coordenadas que permiten visualizar larelación entre unos y otros de forma directa. Este hecho simplifica mucho los algoritmosrespecto a lo que sería una mallado no estructurado. En cualquiera de los casos (estructurado o no estructurado), el mallado deberíacumplir una serie de requisitos genéricos: - Ajuste a las fronteras de la región a estudiar de forma que las condiciones de contorno queden representadas con la mayor exactitud posible. - El mallado debe distribuirse localmente de la forma más regular posible, con variaciones suaves de densidad. La densidad de un mallado se define como el número de puntos por unidad de superficie o volumen. - La mayor densidad del mallado se debe localizar donde se espere que las variaciones espaciales de la solución sean mayores. - El mallado debería ajustarse dinámicamente a las variaciones de las variables en la solución del flujo.Los métodos de generación de mallados utilizados en la práctica se pueden agrupar enlas siguientes categorías (Deconinck, 1994): a) Mallados estructurados. Pueden ser curvilíneos generalizados u ortogonales.Los primeros son más fáciles de construir pero tienen el inconveniente de generar términosadicionales a la hora de resolver las ecuaciones de gobierno. Existen dos grandes grupos detécnicas de construcción de mallados curvilíneos generalizados, la interpolación algebraica y 100. 100Técnicas Numéricas en Mecánica de Fluidosla resolución de ecuaciones diferenciales. En cuanto a los ortogonales, estos pueden serconformes (resolución de ecuación diferencial elíptica) o generados resolviendo una ecuacióndiferencial hiperbólica (técnica de avance en el espacio, es decir, válida para flujos externosdonde una frontera es libre). A priori, se desearía la ortogonalidad del mallado, no sólo por lasimplicidad a la hora de resolver las ecuaciones sino porque cualquier desviación del malladorespecto a la ortogonalidad introduce errores de truncatura en cualquier modelo de solución delas ecuaciones de gobierno. b) Mallados no estructurados. Existen dos técnicas básicas: el avance frontal yla triangularización de Delaney. La primera consiste en la construcción de triángulos otetraedros partiendo de la definición de la frontera del dominio hasta completar el mismo. Latriangularización de Delaney parte de una serie de puntos o nodos y los triángulos o tetraedrosse construyen de forma que dichos puntos son los correspondientes centros de masas. c) Técnicas de multibloque. Combinan mallados estructurados y noestructurados. Se crean bloques cuya estructura general sigue una disposición estructurada demallados elementales que pueden tener una disposición interna en mallados estructurados ono estructurados. No se trata de otro tipo de mallado, sino una mezcla de los dos tipos básicos.5.1.- Clasificación de los mallados basada en la conectividad y estructura dedatos. A continuación, se muestran ejemplos de los distintos tipos de mallados utilizados enla práctica, señalándose las ventajas e inconvenientes de cada uno. Dada la dificultad paramostrar mallados tridimensionales, todos los ejemplos son de aplicaciones bidimensionales. Como grandes grupos de mallados existen los estructurados y los no estructurados, enfunción de las relaciones existentes entre las variables en los distintos puntos. En el tipo de mallados estructurados se subdivide la región a estudiar en rectángulos ocubos (según se trate de problemas bidimensionales o tridimensionales). Se muestra un ejemplo en la figura 72. Sus características más relevantes son lassiguientes:• No es necesario dar la conectividad explícitamente en la formación del malladoya que los nodos vecinos son conocidos.• En general, esto significa que el dominio es tranformado sobre un rectángulo oun cubo subdividido uniformemente en celdas. Por lo tanto, cada nodo esidentificado por los índices (i, j) en el caso de dos dimensiones y con (i, j, k) parael caso de tres dimensiones. En dos dimensiones, los nodos vecinos serían (i-1, j ),(i+1, j), (i, j-1), etc.• Puede ser necesaria una transformación del dominio, tal y como se muestra enla figura 72. 101. Técnicas Numéricas en Mecánica de Fluidos 101Figura 72.- Ejemplo de mallado bidimensional estructurado.Ventajas: - La conectividad ordenada, conduce a resoluciones más simples y, es ideal para computadores vectoriales. - Se requiere de menos memoria para el almacenamiento de las variables. - Los métodos implícitos utilizan la estructura del mallado: ADI, métodos de relajación lineal, relajación en el plano. - La suavización y ortogonalidad pueden ser controlados fácilmente.Desventajas: - El manejo de geometrías complicadas no es flexible. - La adaptibilidad solo es posible añadiendo o moviendo líneas de mallado (2D) o superficies de mallado (3D), lo cual no lo hace flexible. - El movimiento de contornos es difícil de manejar.Como aplicación práctica, se puede recordar el ejemplo de la transmisión de calor enuna placa rectangular, ya estudiado en el apartado anterior. Como ilustración de la mallautilizada, se muestra ésta en la figura 73. 102. 102Técnicas Numéricas en Mecánica de FluidosFigura 73.- Mallado bidimensional estructurado.En el tipo de mallados no estructurados se subdivide la región a estudiar en elementosirregulares bi o tridimensionales. En los siguientes apartados se muestran distintos ejemplosprácticos de este tipo de mallados.5.2.- Métodos de generación de mallados no estructurados. Se muestra un ejemplo en la figura 74 Sus características más relevantes son lassiguientes: • La conectividad tiene que ser expresada explícitamente. Por ejemplo, en la figura 74, se tiene que la celda 10 está constituida por los nodos 7, 14 y 37. • Las celdas pueden ser polígonos cuadriláteros, triángulos, etc. 103. Técnicas Numéricas en Mecánica de Fluidos103Figura 74.- Mallado bidimensional no estructurado.Ventajas: - Flexibilidad para manejar geometrías complicadas, adaptabilidad, movimiento de contornos. - La generación de mallado automática es más fácil, incluso para geometrías de tres dimensiones muy complicadas como el generador de mallado ‘mínima entrada utilizada’.Desventajas: - Requiere mayor almacenamiento de variables y precisa un direccionamiento indirecto. - Las resoluciones son más complicadas en general. - La exactitud es más baja, en general, debido a la falta de suavidad del mallado.Como aplicación práctica, se puede mostrar el mallado para la resolución del flujo enuna bomba centrífuga (figura 75). 104. 104 Técnicas Numéricas en Mecánica de FluidosFigura 75.- Mallado bidimensional no estructurado para una bomba centrífuga.5.3.- Mallados “Multiblock”. Consiste en una combinación de mallados estructurados y no estructurados: en lamacro-escala, el mallado está compuesto de bloques no estructurados (macro-elementos). Concada uno de los “macro-elementos” o “bloque”, el mallado es estructurado. Se muestra unejemplo de geometría difícil de mallar con celdas estructuradas en la figura 76 y la soluciónmultiblock en la figura 77.Figura 76.- Geometría altamente irregular difícil de mallar con celdas estructuradas. 105. Técnicas Numéricas en Mecánica de Fluidos105En este tipo de mallados, las coordenadas de cada nodo llevan una componenteadicional que indica el bloque al que pertenecen. Así, para un punto (i,j), se tendría (i,j,l),indicándose que pertenece al bloque “l” y para un punto en una geometría tridimensional(i,j,k), se pasaría a (i,j,k,l). Figura 77.- Mallado multiblock para la geometría de la figura 76.Se trata de solventar algunas de los puntos débiles del trabajo con los malladosestructurados de bloque único.Ventajas: - Se pueden manejar geometrías más complicadas, lo que aumenta la flexibilidad en general. - Permite maneras fáciles de paralelizar en un multiprocesador. - Ahorra memoria en máquinas de secuencias, permitiendo bloques de Jacobi. (Tratamiento secuencial de cada bloque).Desventajas: - La generación de mallados es difícil: La forma de especificar los contornos entre bloques de forma adecuada exige alto grado de experiencia. - Es todavía menos flexible que los mallados sin estructurar, cuando existe movimiento de contornos. 106. 106Técnicas Numéricas en Mecánica de Fluidos5.4.- Mallados ajustados a los contornos (“Body Fitted Coordinates” o BFC).Respeta los contornos de la geometría que va a ser discretizada, tal y como se muestraen las figuras 78 y 79. Figura 78.- Mallado estructurado que se ajusta al contorno de un álabe.Figura 79.- Mallado no estructurado que se ajusta al contorno de un álabe.5.1.5.- Mallados no ajustados a los contornos. El contorno del dominio no es parte del mallado (ver figura 80). Las características más importantes de este tipo de mallado son las siguientes: • El no ajuste a la geometría del problema conlleva en general a una resolución particular para la gran cantidad de discretizaciones de contornos. • Es muy eficiente para campos internos. • Permite fácil adaptabilidad. 107. Técnicas Numéricas en Mecánica de Fluidos 107 Figura 80.- Mallado estructurado que no se ajusta al contorno superficial (de un álabe). En la práctica se han utilizado muchos tipos de mallados, que resultan de lacombinación de las dos clasificaciones propuestas. Además existen tipos intermedios ysoluciones “imaginativas” que han permitido resolver situaciones geométricas complejas, querequieren la aplicación de técnicas altamente sofisticadas. Dado el enfoque que se pretendeaquí, no se puede profundizar en todo este tipo de técnicas, que corresponderían a un cursomucho más avanzado que el aquí propuesto. 108. 108 Técnicas Numéricas en Mecánica de Fluidos 109. Técnicas Numéricas en Mecánica de Fluidos 1096.- Bibliografía.Existe una multitud de publicaciones relacionadas con los métodos numéricos. Entreellas, destacan el libro escrito por Peyret et al. (1983), los de Hirsch (1988), los de Fletcher(1988), el de Hoffman (1989) y, más recientemente, el editado por Gunzburger (1993), que secentra en métodos apropiados para estudio de flujos incompresibles, el de Anderson (1995) yel de Versteeg et al. (1995), sobre el método de volúmenes finitos. Cada uno de losmencionados presenta sus propias peculiaridades, pero entre todos completan una panorámicabastante amplia de un campo continuamente cambiante. A continuación se detallan los libros, artículos y otras publicaciones que han servidode apoyo para la elaboración de este texto:[1] Abdallah, S., Smith, C.F.,“Three-dimensional solutions for inviscid incompressible flow inturbomachines”, Journal of Turbomachinery. Vol. 112, págs. 391-398. July 1990.[2] Anderson, H.H., “Computational fluid dynamics: the basics with applications”. Mc GrawHill, 1995.[3] Aris, R., “Vectors, tensors and the basic equations of fluid mechanics”, Dover publications.New York, 1962.[4] Arnone, A., Stecco, S.S., “Multigrid calculation of incompressible flows for turbomachineryapplications”, XIV IAHR Congress. ISBN 84-7790-101-5. Págs. c-361 a c-368. Madrid. 1991.[5] Arnone, A., Swanson, R.C., “A Navier-Stokes solver for turbomachinery applications”Journal of Turbomachinery. Vol. 115, págs. 305-313. April 1993.[6] Ballesteros, R., “Métodos numéricos aplicados a la mecánica de fluidos y la transmisión decalor”. Curso de doctorado de la Universidad de Oviedo (publicación interna). 1995.[7] Combes, J.F., “Calcul par elements finis de leculement 3D turbulent dans une pompecentrifuge”, AGARD Conference Proceedings 510. CFD Thechniques for PropulsionApplications. San Antonio, USA. 1991.[8] Deconinck, H., “Overview of grid generation methods”, Lecture Series del von KarmanInstitute for Fluid Dynamics. Bélgica, 1994.[9] Degrez, G., “Introduction to computational fluid dynamics”, Lecture Series del von KarmanInstitute for Fluid Dynamics. Bélgica, 1994.[10] Fletcher, C.A.J., “Computational techniques for fluid dynamics” (vols. I y II). Springer-Verlag, USA. 1988. 110. 110Técnicas Numéricas en Mecánica de Fluidos[11] Fox, R.W., McDonald, A.T. “Introducción a la mecánica de fluidos”. McGraw-Hill1,Méjico 1995 (cuarta edición).[12] Gresho, P.M., “Incompressible fluid dynamics: some fundamental formulation issues”,Annual Review of Fluid Mechanics. Vol. 23, págs. 413 a 453. 1991.[13] González, J., Santolaria, C., Ballesteros, R., Robles del Peso, A., “Modelización numéricade flujos incompresibles no estacionarios mediante un método de paso fraccional”, Anales deIngeniería Mecánica. Año 11, vol. I, págs. 35-42, 1997.[14] González, J., Santolaria, C., Ballesteros, R., “Desarrollo de un algoritmo para la resoluciónde flujos no estacionarios bidimensionales usando un método de paso fraccional”, Anales deIngeniería Mecánica. Año 12, vol. I, págs. 553-558, 1998.[15] González, J., “Modelización numérica del flujo no estacionario en una bomba centrífuga.Efectos dinámicos de la interacción entre rodete y voluta”, Tesis Doctoral, Universidad deOviedo, 2000 (publicada en CD ROM).[16] Gunzburger, M.D., Nicolaides, R.A., “Incompressible computational fluid dynamics.Trends and advances”, Cambridge University Press. 1993.[17] Hirsch, C., “Numerical computation of internal and external flows. Vols. I y II:fundamentals of numerical discretization & computational methods for inviscid and viscousflows”, John Wiley & Sons. Bélgica, 1988.[18] Hirsch, C., “State of the art of computational fluid dynamics in industry”, IndustrialComputational Fluid Dynamics. Lecture Series 1995-03. Von Karman Institute, Bélgica. 1995.[19] Hoffman, K.A., “Computational fluid dynamics for engineers”, Engineering EducationSystem. Austin, USA. 1989.[20] Jameson, A., Schmidt, W., Turkel, E., “Numerical solutions of the Euler equations by finitevolume methods using Runge - Kutta time - stepping”, AIAA Paper No. 81-1259. 1981.[21] Jameson, A., “The role of CFD in preliminary aerospace design”, ASME FEDSM2003-45812, Honolulu, (2003).[22] Kwak, D., Chakravarthy, S.R., “A three-dimensional incompressible Navier-Stokes flowsolver using primitive variables”, AIAA Journal. Vol. 24, No. 3, pags. 390-396. March 1986.[23] Kwak, D., Kiris, C., Wiltberger, N., Rogers, S, Rosenfeld, M., "Numerical methods forsimulating unsteady incompressible flows", Lecture Notes in Phisics, Vol. 414, págs. 448-452.Numerical Methods in F.D., 13th Conference. Rome, 1992.[24] Lakshminarayana, B., “An assesment of computational fluid dynamic techniques in theanalysis and design of turbomachinery – The 1990 Freeman Scholar Lecture”, Journal of FluidEngineering, vol. 113, pp. 315-352. Septiembre 1991.[25] Launder B.E., Spalding D.B., “Mathematical models of turbulence”, Academic Press. NewYork, 1972. 111. Técnicas Numéricas en Mecánica de Fluidos 111[26] P.L. Lions, “Mathematical topics in fluid mechanics. Incompressible models”. ClarendonPress, Oxford, 1996.[27] Moin, P., Kim, J., “Tackling turbulence with supercomputers”, Scientific American, 276.1997.[28] Niederdrenk, P., “A short introduction to numerical grid generation”, Course Note (IB 221-87 A 14) at DFVLR-AVA, Göttingen. 1987.[29] Pan, D, Chackravarthy, S.R., “Unified formulation for incompressible flow”, AIAA Paper,No. 89-01222. 27th Aerospace Science Meeting. Reno, Nevada. 1989.[30] Panton, R.L., “Incompressible flow”, John Wiley & Sons, USA, 2nd Edition. 1996.[31] Patankar, T.C., “Numerical heat transfer and fluid flow”, Mc. Graw-Hill. 1980.[32] Pentaris, A., Tsangaris, S., “A proyection methodology for the simulation of unsteadyincompressible viscous flows using the approximate factorization technique”, NATO-AGARDCP-578. Progress and Challenges in CFD Methods and Algorithms, págs. 35-1 a 35-14. Abril,1996.[33] Peyret, R., Taylor, T.D., “Computational methods for fluid flow”, Springer-Verlag. NewYork, 1983.[34] Quarteroni, A., “Domain decomposition methods: algebraic aaspects and application tocompressible and incompressible flows”, 27th Computational Fluid Dynamics. Lecture Series.Von Karman Institute. Bélgica. 1996.[35] Rasmussen, E.B., “A finite diference scheme for three-dimensional modelling of fluiddynamics”, XIV IAHR Congress. ISBN 84-7790-101-5. Págs. c -341 a c - 348. Madrid. 1991.[36] Shaw, C.T., “Using computational fluid dynamics”, Prentice Hall. 1992.[37] Strazisar, A.J., “The changing roles of experimental and computational fluid mechanics”,Global Gas Turbines News, IGTI, pp. 16-18. Agosto 1994.[38] Streeter, V.L., Wylie, E.D., “Mecánica de los Fluidos”, Mc. Graw-Hill, Méjico. 1987.[39] Thomson, J.F., Warsi, Z.U.A., Mastin, C.W., “Numerical grid generation: foundations andapplications”, Ed. Joe F. Thomson. North-Holland, 1982.[40] Thomson, J.F., “Grid generation techniques in computational fluid dynamics”, AIAAJournal, Vol. 22, No. 11, pp. 1505-1523. 1984.[41] Versteeg, H.K., Malalasekera, W., “An introduction to computational fluid dynamics”,Longman Scientific&Technical. 1995.[42] White, F.M. “Mecánica de fluidos”. McGraw-Hill, Madrid 1983. 112. 112Técnicas Numéricas en Mecánica de Fluidos[43] White, F.M., “Viscous fluid flow”. Mc. Graw-Hill Editions. 1991.[44] Wilcox, D.C., “Turbulence modeling for CFD”, DCW Industries. 1993. 113. Técnicas Numéricas en Mecánica de Fluidos113Anexo I. El metodo Von Neumann paraanalisis de estabilidad.Los análisis de estabilidad para problemas lineales con coeficientes constantes sonbien conocidos si las condiciones de contorno pueden ser despreciadas. Este es el caso para undominio infinito, o para condiciones periódicas en un dominio finito. En este ultimo caso seconsidera que la longitud L en el eje x se repite periódicamente, por tanto, todas lascantidades, la solución, así como también los errores, pueden ser desarrollados en una serie deFourier finita sobre el dominio 2L. Este desarrollo forma las bases para el análisis deestabilidad de Von Neumann. Este método es el más utilizado para los análisis de estabilidad en los que intervienenecuaciones lineales con coeficientes constantes y condiciones periódicas de contorno. Sinembargo, siempre que se tenga que tratar con coeficientes no constantes, y términos nolineales en las ecuaciones básicas, la información sobre la estabilidad se vuelve muy limitada.Consideraciones iniciales. Considerando uno de los modelos representativos de la convección, se considera lasiguiente ecuación hiperbólica: ∂u∂u+a= 0(AI.1) ∂t∂xdonde “u” función de ( x,t ); “a” la velocidad de convección, o la velocidad de onda, deacuerdo a la interpretación dada a la expresión AI.1 y se supone que a > 0. Se suele utilizaruna notación más corta de la misma, donde la derivada esta indicada como subíndice: ut + a ux = 0 (AI.2) Se consideran las siguientes condiciones iniciales: t=0u ( x, t ) = f ( x ) 0≤x≤L (AI.3) x=0u (0, t ) = g ( x )t≥0 Para aplicar un método de elementos finitos a la ecuación, se podría establecer uncentro en “i”. Haciendo la subdivisión del espacio del dominio en celdas de longitud ∆x , laformula de discretización de u x para el punto “i” del mallado es la siguiente:a (u t ) i = −(u i +1 − u i −1 )(AI.4)2 ∆x El lado izquierdo representa la derivada respecto al tiempo en el punto “i”. 114. 114Técnicas Numéricas en Mecánica de FluidosEl siguiente paso es la discretización del tiempo. Esto implica el remplazo de (u t ) ipor la forma discreta, teniendo en cuenta el nivel del tiempo para el cual el lado derecho del laecuación será evaluado. Un esquema simple sería la evaluación del lado derecho en intervalosde tiempo “n”. Tal método es conocido como el método de Euler para la integración respectoal tiempo de una ecuación diferencial ordinaria, que es un esquema explícito.u in +1 − u ina =−(u in+1 − u in−1 ) (AI.5)∆t2 ∆x Evaluando el lado derecho de la ecuación en el nivel (n+1) conduce al esquemaimplícito, conocido como aguas abajo, o método de Euler implícito u in +1 − u ina=−(u in++11 − u in−+11 )(AI.6) ∆t2 ∆xdonde tres valores aparecen simultáneamente en el nivel (n+1) De las definiciones de orden de precisión de las fórmulas de diferencias finitas seespera que las ecuaciones AI.5 y AI.6 sean de primer orden en tiempo, y de segundo orden enespacio en los puntos (“i”) y en el nivel de tiempo “n”. Por otro lado, mediante unadiferenciación en espacio aguas abajo, con una aproximación de primer orden, se llega a laforma semidiscreta:a(u t ) i = −(u i − u i −1 ) (AI.7) ∆xCon una diferenciación hacia delante en tiempo se obtendría el siguiente esquemaexplícito, conocido como el esquema de primer orden aguas arriba (“upwind”), es decir: u in +1 − u in a=−(u in − u in−1 )(AI.8) ∆t∆x La correspondiente versión implícita, evaluando el lado derecho (n+1) es:u in +1 − u in a =−(u in +1 − u in−+11 )(AI.9)∆t∆xa ∆tdonde se define el cociente σ =, que constituye el denominado habitualmente como ∆x“número de Courant”. 115. Técnicas Numéricas en Mecánica de Fluidos115Descomposición del error (Fourier)n Sea u i la solución exacta de una ecuación, y u in la solución numérica aproximada. Ladiferencia podría ser debida a los errores de redondeo introducidos por el algoritmo decálculo, o bien a los errores en los datos iniciales. De esta forma, se debe cumplir: nu in = u i + ε in o bien, lim ε in ≤ K(AI.10) n →∞donde ε in indica el error respecto al tiempo en el punto “i” del mallado. K es independiente denn. Cualquier esquema numérico lineal es satisfecho exactamente por u i , y por tanto loserrores ε in son también soluciones para la ecuación. Con el fin de presentar los fundamentos del método, primero se hará referencia a losejemplos previos. Considerando la ecuación AI.5, e insertando en ella el resultado de AI.10 setiene: n +1ui− u in ε in +1 aa nn+=− (ε in+1 − ε in ) − (u i +1 − u i −1 ) (AI.11) ∆t ∆t2∆x2∆xn Como u i satisface la ecuación AI.5, se obtiene la ecuación para los errores ε in :ε in +1 − ε in a =− (ε in+1 − ε in−1 )(AI.12)∆t2∆xque es idéntica al esquema básico. Así pues, los errores ε in varían a través del tiempo de lamisma forma que la solución numérica u in .Los errores no deberían crecer indefinidamente de un intervalo de tiempo a otro. Unanálisis basado en el desarrollo del tiempo en si mismo en lugar de la conducta del error,establece que cualquier componente de la solución inicial no debería ser amplificado sinlimite. Expresada esta condición en forma matemática se tiene una matriz en la cual todos lasincógnitas de cada punto del mallado, en un tiempo dado, están agrupadas en un vector “U”,definido como sigue un tiempo n ∆ t . u1  n ...  n u i −1 U = n n(AI.13)u  i  u in+1  ...  El esquema se puede escribir con un operador como función de Un según la expresión: 116. 116Técnicas Numéricas en Mecánica de Fluidos n +1 nU =CU(AI.14)donde el operador “C” es función del tiempo ∆ t , y del tamaño del mallado ∆ x . Si se considera al operador “C”, lineal, si ε n designa el vector columna de los erroresen un nivel de tiempo “n”:   εni−1  εn=  ε n  (AI.15) i ε ni+1  nLa ecuación AI.10 se puede escribir con U indicando la solución exacta,Un = U n + en(AI.16)Insertando esta ecuación en el esquema básico se conduce a la expresión: n +1 n U+ e n +1 = C U + C e n (AI.17)o bien:e n +1 = Ce n(AI.18) npor definición de U como una solución de la expresión:n +1nU =CU(AI.19) Así, la evolución en el tiempo del error esta determinado por el mismo operador Ccomo solución numérica del problema. Si las condiciones de contorno están consideradascomo periódicas, el error ε in puede ser descompuesto en series de Fourier en el espacio paracada nivel de tiempo “n”. Como el dominio es de una longitud finita, se tendrá unarepresentación de Fourier sobre un número finito de oscilaciones. En el dominio L la serie deFourier representa la región (0, L) y la negativa (-L, 0),la frecuencia fundamental correspondea la máxima longitud de onda λmax = 2L. La longitud de onda asociada k = 2π/λ alcanza unmínimo valor en kmin = π/L. Además, el máximo valor de kmax en el intervalo (-L, L) estáasociado con la menor longitud de onda con un espaciamiento ∆x. Esta menor longitud deonda es igual a 2∆x (ver figura AI.1) y por tanto, kmax = π/∆x.Así, con un mallado de índice “i”, de 0 a N, con xi = i. ∆x, se tendrá: 117. Técnicas Numéricas en Mecánica de Fluidos 117∆x = L/N (AI.20)Todas las oscilaciones representadas en un tiempo finito están dadas por:π π k j = j k min = j=j con j = 0,1,2,...,N (AI.21)LN ∆xcon el máximo valor de j asociado a la máxima frecuencia. Así, con kmax = π/∆x el mayorvalor de “j” es igual al número de intervalos N del mallado (ver figura AI.1). Figura AI.1.- Representación del error en el intervalo (-L, L).Cualquier función finita del mallado se puede descomponer en series de Fouriercomo:NN ∑ E nj e ∑E I k j i ∆xε in ==n j eI i j π/ N (AI.22) j= − N j= − Ndonde I = − 1 y E in es la amplitud de la oscilación “j-ésima”. Por otro lado, j = 0 representa una función constante. El producto kj ∆x es a menudorepresentado como un ángulo de fase. πΦ = k j ·∆x = j(AI.23) N 118. 118Técnicas Numéricas en Mecánica de Fluidos Y cubre el dominio (-π, π) en intervalos π/N.Factor de amplificación.El tiempo de evolución de una oscilación simple E n e IiΦ es determinado de la misma formajque u in . Así utilizando estos términos en la ecuación AI.12, se obtiene:( E n +1 + E n )a e Ii Φ+ ( E n e I ( i +1 ) Φ − E n e I ( i −1 ) Φ ) = 0 (AI.24) ∆t2∆x Dividiendo entre e Ii Φ , se llega a: n +1σ n IΦ(E+ En)+ E (e − e − I Φ ) e Ii Φ = 0 (AI.25) 2donde el parámetro σ esa ∆tσ= (AI.26)∆x La condición de estabilidad dada por AI.10 será satisfecha si la amplitud de cualquiererror En no crece en el tiempo, es decir, si la relación:E n +1 G ≡ ≤1 para todo Φ(AI.27)En La cantidad G, definida porE n +1 G=(AI.28) Enes el denominado factor de amplificación, y es función del factor del intervalo de tiempo ∆t, afrecuencia y el valor de ∆x. De la ecuación AI.25 se puede llegar a:σG −1+ 2 I sin Φ = 0(AI.29)2o bien, una expresión equivalente sería:G = 1 − I σ sin Φ(AI.30) La condición de estabilidad requiere que G sea menor o igual a uno: 119. Técnicas Numéricas en Mecánica de Fluidos1192G = 1+ σ 2 sin 2 Φ(AI.31)y esto claramente no es satisfecho.Ejemplo del esquema dado por la discretización AI.8: estabilidad condicional.Insertando la oscilación simple E n e IiΦ en el esquema dado por la ecuación AI.8, sejobtiene:n +1 (E− E n ) + σ E n ( e Ii Φ − e − I ( i − 1 ) Φ ) = 0 (AI.32)O bien, después de dividir por E n e IiΦ j G = 1 − σ + σe − IΦ(AI.33)es decir: G = 1 − 2 σ sin 2 Φ / 2 − Iσ sin Φ (AI.34)Con el fin de analizar la estabilidad de la ecuación AI.8, que es la región donde elmodulo del factor de amplificación G es menor a uno, una conveniente aproximación es unarepresentación de G en el plano complejo. Escribiendo ξ y η, respectivamente para las partesreal e imaginaria de G se tendrá:2ξ = 1 − 2 σ sin Φ / 2 = (1 − σ ) + σ cos Φ(AI.35)η = − σ sin Φlo cual puede ser considerado como una ecuación paramétrica para G con Φ como parámetro.En el plano complejo de G la condición de estabilidad dada por AI.27 fija que la curvarepresentativa de G para todos los valores de Φ = k·∆x debería permanecer dentro del círculo(ver la Figura AI.2). El esquema es estable para:0 < σ ≤1(AI.36) 120. 120Técnicas Numéricas en Mecánica de Fluidos Figura AI.2.- Región de estabilidad para el esquema dado por la expresión AI.8.Así la ecuación AI.8 es condicionalmente estable, y la condición AI.36 es conocidacomo la Courant-Friedrich-Lewy o condición CFL. El parametro σ es llamado el número deCourant.Comentarios sobre la condición CFL.Esta condición de estabilidad de la mayoría de los esquemas explícitos para ondas yecuaciones de la convección expresa que la distancia cubierta durante el intervalo de tiempo“∆t”, por los disturbios de propagación a velocidad “a”, deben ser más bajos que la distanciamínima entre dos puntos del mallado. Teniendo en cuenta la figura AI.3, la línea PQ es lacaracteristica dx/dt = a (a través de P) y define el dominio de dependencia de la ecuacióndiferencial en P. Por otro lado la ecuación define el dominio de dependencia de P entre PAC.La condición σ ≤ 1 expresa que ∆t/∆x debe ser elegido de tal forma que el dominio dedependencia de la ecuación diferencial debería estar contenido en el dominio de dependenciade la ecuación discretizada. Es decir, el esquema numérico u in +1 en el punto i debe ser capazde incluir toda la información física que influye en ese punto. 121. Técnicas Numéricas en Mecánica de Fluidos 121 Figura AI.3.- Interpretación geométrica de la condición CFL, σ ≤ 1 .Ejemplo del esquema dado por AI.6: condición de inestabilidad.El esquema implícito de Euler, con una discretización central del espacio de laecuación de convección ofrece una tercera situación con respecto a las condiciones deestabilidad. Con el mismo análisis de estabilidad que AI.6, la amplitud del error En+1 sevuelve, después de la introducción de una oscilación de la forma En eΙiΦ σ e Ii Φ ( E n +1 + E n ) + E n +1( e I Φ − e − I Φ ) e Ii Φ = 0 (AI.37) 2o bien:σ G −1+G ( e IΦ − e − IΦ ) = 0(AI.38)2y se tiene1 G = (AI.39) 1 + I σ sin ΦEl modulo G siempre es menor que uno para todos los valores de σ 2 1 G = G .G * =2 2 (AI.40)1+ σ sin Φ Por tanto la ecuación implícita es incondicionalmente estable. Así, esta visto que elsistema puede tener estabilidad condicional, estabilidad incondicional, o inestabilidadincondicional. El método Von Neumann ofrece una forma fácil y simple de hallar laspropiedades de estabilidad de los esquemas lineales con coeficientes constantes cuando lascondiciones de frontera son asumidas como periódicas. 122. 122Técnicas Numéricas en Mecánica de FluidosFormulación general del metodo Von Neumann: sistema de ecuaciones. Se muestra aquí la formulación en forma de matriz y operador. Se considera que elesquema numérico es obtenido en dos pasos: un espacio de discretización, seguido por eltiempo de integración.1) Cuando se aplica el espacio de discretización (método de diferencias finitas) el espacio diferencial del operador es aproximado por un espacio discretizado del operador S, conduciendo al método de formulación de línea para valores discretos u in = u ( x , n∆t ) … donde … xi es el punto "i" del mallado. du i= S ui + qi (AI.41)dtdonde qi: contiene términos fuentes eventuales y la contribución desde las condiciones defrontera. La representación matricial es descrita con el vector Un del sistema dado por AI.15.dU =SU+q(AI.42)dt2) Cuando se aplica el esquema de tiempo de integración, correspondiendo a un nivel dos del esquema conectando niveles de tiempo (n+1) y n, el esquema numérico asociado con lan +1 ecuación diferencial genérica, es decir: U= C U n , se tiene: nu n +1 = C u i + q i(AI.43)o en forma matricial:U n +1 = C U n + Q(AI.44)donde C puede ser considerado como un operador de discretización del esquema. Para el nivel dos del esquema implícito de la forma: B1 U n +1 = B 0 U n el operador C esdefinido por C = B1−1B 0 . 123. Técnicas Numéricas en Mecánica de Fluidos123Anexo II. Ecuacionesen sistemasgeneralizados. Partiendo de una ecuación diferencial en la forma∂ q∂ F ∂ G + +=0 ∂ t ∂ x ∂ y(AII.1)que se puede escribir en notación simplificada: q t + Fx + G y = 0(AII.2) Si se realiza el cambio de variables independiente del tiempo, pasando a un nuevosistema de coordenadas ξ, η: ( ξ, η ) = f ( x, y ) ≠ f ( t )(AII.3) Para modificar las derivadas parciales se debe hacer siguiendo la regla de la cadena, esdecir: ∂ ∂∂ = ξx+ ηx ∂ x ∂ ξ∂ η(AII.4) ∂ ∂∂ = ξy+ ηy ∂ y ∂ ξ∂ η Según este cambio de variables, el Jacobiano de la transformación resulta ser:∂ ξ∂ ξ∂ x∂ y ξx ξy J== = ξ x ηy − ξ y ηx ∂ η ∂ η ηx ηy∂ x∂ y(AII.5) Y, entonces, la ecuación diferencial de partida se transforma en: 124. 124 Técnicas Numéricas en Mecánica de Fluidos∂ q∂ F ∂ F ∂ G ∂ G+ ξx + ηx + ξy+ ηy =0∂ t∂ ξ ∂ η ∂ ξ ∂ η(AII.6) Ecuación que en notación simplificada resulta: q t + ξ x Fξ + ξ y G ξ + η x Fη + η y G η = 0(AII.7) Dividiendo por el Jacobiano, que es independiente de la variable tiempo, se obtiene:∂  q  ξ x Fξ ξ y G ξ η x Fη η y G η  + + ++=0∂ tJ JJ JJ(AII.8) A partir de operaciones para transformar las parciales con respecto a (x,y) a parcialesrespecto a las nuevas variables (ξ, η), se obtiene finalmente: ∂ q* ** + Fξ + G η = 0 ∂ t(AII.9)donde las variables con asterisco representan: ξx F + ξyG ηx F + ηy GF* =, G* = J J (AII.10)y q* representa la variable q dividida por el Jacobiano de la transformación.Ecuaciones para flujo incompresible ideal bidimensional. Ecuación de continuidad: ∂ u∂ v+ =0 ∂ x∂ y (AII.11) Ecuación de cantidad de movimiento:∂ u ∂ p∂ u2∂ ( u v)++ +=0∂ t ∂ x∂ x ∂ y(AII.12)∂ v ∂ p∂ ( u v)∂ v2++ +=0∂ t ∂ y∂ x ∂ y 125. Técnicas Numéricas en Mecánica de Fluidos125 Adimensionalizando según las variables características del problema de descarga deldepósito, se tiene: xx* = L yy* = Lu(AII.13)u* = 2gHvv* = 2gHLas ecuaciones no sufren variación y serían las mismas pero con las variablesadimensionales.Cambio de sistema de coordenadas (plano físico-plano computacional). Se define el cambio de coordenadas siguiente:ξ = x *(AII.14) 2 y*η =* * * h 1 + (h 2 − h 1 ) x * Es decir, prescindiendo de los asteriscos, se tendrán las siguientes métricas:∂ ξ=1∂ x∂ ξ=0∂ y (AII.15)∂ η− 2 (h 2 − h 1 ) y (h 2 − h 1 ) η= =∂ x [h 1+ (h 2 − h 1 ) x] 2[h 1+ (h 2 − h 1 ) ξ ]∂ η−2−2==∂ y [h 1 + (h 2 − h 1 ) x] [ h 1 + (h 2 − h 1 ) ξ] Y el Jacobiano, por tanto, resulta ser: 126. 126 Técnicas Numéricas en Mecánica de Fluidos ξx ξy 2 J== ξ x ηy − ξ y ηx = ηx ηy [ h 1 + (h 2 − h 1 ) ξ ](AII.16) La ecuación de continuidad quedará en este caso, según la ecuación (AII.9): ∂ u 2 (h 2 − h 1 ) y ∂u2 ∂v −+ =0 ∂ ξ (h 1 + (h 2 − h 1 ) x )2∂ η (h 1 + (h 2 − h 1 ) x )2∂ η (AII.17) Es decir: ∂ u∂u ∂v +A +B =0(AII.18) ∂ ξ∂ η∂ η Con las constantes A y B iguales a: 2 (h 2 − h 1 ) y A=−(h 1 + (h 2 − h 1 ) x )2 (AII.19) B=J Las ecuaciones de cantidad de movimiento (AII.12) se transforman según la expresióngeneral (AII.9), es decir:∂ ξ∂ ξ ∂ η∂ η( p + u2 ) + ( u v)  ( p + u 2) + ( u v) ∂ u ∂ ∂ x∂ y∂∂ x∂ y∂ + + =0t ∂ ξJ ∂ η  J  ∂ ξ ∂ ξ∂ η ∂ η  ( u v) + ( p + v 2) ( u v) +( p + v 2) ∂ v ∂  ∂ x∂ y ∂∂ x ∂ y+  +   =0∂ t ∂ ξ J ∂ η J      (AII.20) Operando, para el caso particular del cambio de variables definido, se llega a lasexpresiones: 127. Técnicas Numéricas en Mecánica de Fluidos 127∂ u∂  p + u2∂  − 2 (h 2 − h 1 ) y p + u2 2uv+  + ∂ η  h + (h − h ) x +=0∂ t ∂ ξ J   1(2 1 )2 J ( h 1 + (h 2 − h 1 ) x J )∂ v∂  u v  ∂  − 2 (h 2 − h 1 ) yuv 2( p + v 2) ++  +  =0∂( t ∂ ξ  J  ∂ η  h 1 + (h 2 − h 1 ) x)2 J ( h 1 + (h 2 − h 1 ) xJ) (AII.21) Ecuaciones que se simplifican bastante si se sustituye el valor del Jacobiano según elcálculo realizado en la ecuación (AII.16). Así, se llega a: ∂ u∂ p + u2 ∂  −η ( h 2 − h1 )  +  J+  ∂ η  ( )p + u2 + u v = 0 ∂ t ∂ ξ  2   (AII.22)  ∂ v∂  u v  ∂ − η ( h 2 − h 1 )  ∂+ +( u v) + ( p + v2 )  = 0 t ∂ ξ  J  ∂ η 2 Si ahora se usa la ecuación de continuidad (AII.17), se puede simplificar esta últimaexpresión, obteniéndose:∂ u∂  p   u  ∂ J η ( h 2 − h1 ) ∂ p 2 ( h 2 − h1 )∂ u∂v +  −   − −p + u2 + v −u(=0 )∂ t ∂ ξJ J ∂ ξ2 ∂ η 2 ∂ η∂ η∂ v u∂ v ∂ 1  η ( h 2 − h 1 ) ∂v( h 2 − h1 )∂p ∂v ++ ( u v)   − u −u v+ +v =0∂ t J∂ ξ∂ ξ J2∂ η 2 ∂ η∂ η(AII.23) 128. 128 Técnicas Numéricas en Mecánica de Fluidos 129. Técnicas Numéricas en Mecánica de Fluidos 129Anexo III. Glosario de términos empleadosen técnicas numéricas.Las técnicas numéricas aplicadas a la Mecánica de Fluidos se prestan en la mayoría delos textos a la aparición de multitud de acrónimos o siglas, tanto de distintas técnicas utilizadas,como de instituciones que se dedican a la investigación en este campo. Esto puede sorprenderinicialmente e incluso crear gran confusión. Este anexo surgió con el fin de poder servir de ayudaen la interpretación de alguno de los mencionados términos. No se pretende aclarar todos losmatices de dichos términos (algunos de los cuales darían para escribir libros enteros) sino que setrata de indicar el brevemente su significado y su contexto de aplicación.Se incluye a continuación un glosario de términos empleados habitualmente en textosy artículos relacionados con las técnicas numéricas. Algunos de los términos estánrelacionados con la aplicación de dichas técnicas numéricas al campo más específico de lasturbomáquinas dentro de la Mecánica de Fluidos y algunas son más generales. También seincluyen los acrónimos de algunas instituciones u organismos que aparecen comúnmente enlos textos.La organización de la información se hace dividiendo la lista completa en tres subgrupos,uno con los acrónimos de técnicas, algoritmos o estrategias numéricas, otro con las siglas deinstituciones u organismos y un último con los acrónimos no clasificables en los dos gruposanteriores. Dentro de cada grupo, se ordenan las palabras por orden alfabético de sus iniciales yse explica su significado siguiendo el esquema tipo que se indica a continuación: ACRÓNIMO EN MAYÚSCULAS (significado en inglés de las iniciales). Breveexplicación del origen o significado del término.AIII.1.- Acrónimos usuales en los textos de Mecánica de FluidosComputacional.A AC (Aditive Corrections). Correcciones en los valores de las variables en técnicas de multimallados. De forma más general se usan para indicar corriente alterna. ADI (Alternating Direction Implicit). Método de resolución de ecuaciones lineales. Como indica su nombre se trata de un método implícito según direcciones alternativas. También se aplica este término para la solución de distintos algoritmos que siguen el mencinado procedimiento. AF (Approximate Factorization Scheme). Técnica numérica de descomposición de derivadas por términos equivalentes hasta un determinado orden. AFLBI (Approximate Factorization Linear Block Implicit Scheme). Esquema implícito basado en la factorización aproximada (AF) de las ecuaciones. AMG (Algebraic Multi Grid). Estrategia multimallado algebraica. Constituye uno de los métodos más rápidos para solución de sistemas lineales de ecuaciones. 130. 130Técnicas Numéricas en Mecánica de Fluidos ARSM (Algebraic Reynolds Stresses Model). Modelo de turbulencia algebraico de lastensiones de Reynolds.B BL (Baldwin-Lomax). Modelo de turbulencia algebraico. BC (Boundary conditions). Condiciones de contorno de un determinado problema. BFC (Body Fitted Coordinates). Coordenadas que siguen líneas de los contornos.C CAD (Computer Aided Design). Diseño asistido por ordenador. Conjunto de técnicas de dibujo en ingeniería apoyadas en el uso de ordenadores. CAE (Computer Aided Engineering). Ingeniería asistida por ordenador. Similar al término anterior pero más amplio. CAM (Computer Aided Manufacturing). Fabricación asistida por ordenador. Término también utilizado en contextos muy similares a los dos anteriores. CFD (Computational Fluid Dynamics). Técnicas numéricas aplicadas a la Mecánica de Fluidos. CGC (Coarse Grid Correction). Modificación de la solución en técnicas multimallado a partir del valor de las variables en la malla menos fina. CGM (Congugate Gradients Method). Método de los gradientes conjugados. CPPE (Consistent Pressure Poisson Equation). Método consistente utilizado en la ecuación de Poisson para la presión. El método de la ecuación de Poisson para la presión se explica en el capítulo cuatro de esta tesis. CPU (Central Proccess Unit). Unidad de procesado de datos de un ordenador. Constituye el núcleo básico de cualquier sistema informático. CUPW (Convected Upwind Scheme). Algoritmo de solución de las ecuaciones aguas arriba. Desarrollado por Lacor en 1986.D DAES (Differential Algebraic Equations). Ecuaciones diferenciales resueltas por métodosalgebraicos. DO (Discrete Ordinates Radiation Model). Modelo de radiación discreto. DOC (Direct Operating Costs). Costes de operación directos. Término usado habitualmenteen aeronáutica. DRW (Discrete Random Walk). Movimientos discretos y aleatorios. DS (Direct Simulation). Método de simulación directa de las ecuaciones de gobierno. DTRM (Discrete Transfer Radiation Model). Modelo de radiación de la transferenciadiscreta.F FAS (Full Aproximated Storage). Método de almacenamiento de datos. FAVOR (Fractional Area/Volume Representation). Algoritmo que utiliza una mallacartesiana como base para el dominio de la solución, calculando fracciones devolumen o área para cada cara de la discretización. FCBR (Fully Coupled Blade Row). Tipo de cálculo del flujo en una turbomáquina queconsidera los efectos de varias coronas de álabes. 131. Técnicas Numéricas en Mecánica de Fluidos 131FCT (Flux Corrected Transport). Se trata del método precesor del TVD y fue introducido por Boris y Brook en 1973.FDM (Finite Difference Method). El método de elementos finitos consiste en una discretización numérica de las ecuaciones de gobierno basada en la aproximación de derivadas por diferencias entre las variables.FEM (Finite Element Method). El método de elementos finitos consiste en una discretización numérica de las ecuaciones de gobierno basada en la utilización de funciones de interpolación para las distintas variables.FVM (Finite Volume Method). El método de elementos finitos consiste en una de discretización numérica de las ecuaciones de gobierno basada en la propiedad de conservación de dichas ecuaciones.GGCGTVD (Generalized Conjugate Gradient Total Variation Diminishing). Técnica numérica derivada de los métodos TVD.GGDH (Generalized Gradient Diffusion Hypothesis). Hipótesis del gradient generalizado de difusión. Introducida por Daly y Harlow en 1970.HHOT (Higher Order Terms). Términos de orden superior.HP (High Pressure Stages). Escalones (Stages) de turbomáquinas de alta relación depresiones.IIBVP (Initial Boundary Value Problem). Tipo de problema cuya solución vienecondicionada por una determinada situación inicial o valor inicial de las variables.IC (Initial Conditions). Condiciones iniciales de un problema que evoluciona en el tiempo.IGV (Inlet Guide Vanes). Corona de álabes a la entrada de una turbomáquina.LLED (Local Extremum Diminsishing). Técnica de solución de las ecuaciones desarrollada por Jameson en 1993. Se trata de un esquema no oscilatorio de segundo orden de precisión.LES (Large Eddy Simulation). Aproximación a las ecuaciones de Navier-Stokes usando las hipótesis de intercambio energético entre los vórtices de gran y pequeño tamaño. Discutido más en profundidad en el capítulo cuatro del texto.LHS (Left Hand Side). Lado izquierdo de una igualdad matemática.LP (Low Pressure Stages). Escalones de baja presión en turbomáquinas.LSOR (Line Succesive Overrelaxation). Técnica de resolución de sistemas de ecuaciones lineales basada en la utilización de un parámetro de sobrerrelajación.LU (Lower and Upper Diagonal). Método de resolución de sistemas de ecuaciones lineales basado en los elementos por encima y por debajo de la diagonal.LU-SGS (Lower-Upper Symmetric-Gauss-Seidel Scheme). Similar al anterior pero utilizando una estrategia Gauss-Seidel simétrica.MMAC (Marker and Cell). Modelos de resolución de las ecuaciones diferenciales basados en la utilización de mallados escalonados. 132. 132 Técnicas Numéricas en Mecánica de Fluidos MARS (Monotone Advection and Reconstruction Scheme) Esquema de segundo orden conmenor sensibilidad o precisión con la estructura de la malla y el “skewness”. MHD (Magneto Hidrodynamics). Magnetohidrodinámica de Fluidos. MRF (Multiple Reference Frame). Métodos de resolución de las ecuaciones que se basan enla utilización de varios sistemas de referencia. MUSCL (Monotone Upstream-Centered Schemes for Conservation Laws). Leyes deconservación resueltas según algoritmos monótonos aguas arriba.N NBC (Natural Boundary Conditions). Condiciones de contorno naturales. NS (Navier-Stokes). Ecuaciones de gobierno en Mecánica de Fluidos.O OBC (Outflow Boundary Conditions). Condiciones de contorno de cualquier tipo, excepto del tipo Dirichlet. También a veces se usa para hablar de Open Boundary Conditions, con el mismo significado. ODE (Ordinary Differential Equations). Ecuaciones diferenciales simples, en contraposición con las PDE.P PDE (Partial Differential Equations). Ecuaciones diferenciales en derivadas parciales. PISO (Pressure Implicit with Spliting of Operators). Método de solución de las ecuacionespara flujo incompresible. PNS (Parabolized Navier-Stokes). Modelo de las ecuaciones de Navier-Stokes“parabolizadas”. Modelo intermedio entre el TSL y las ecuaciones reales. PPE (Pressure Poisson Equation). Método de resolución de las ecuaciones de gobierno paraflujo incompresible consistente en operar haciendo la divergencia de lasecuaciones de cantidad de movimiento y luego aplicar la condición deincompresibilidad. PPNS (Partially Parabolized Navier Stokes). Idéntico al PNS pero haciendo la“parabolización” de forma parcial. PSOR (Point Succesive Overrelaxation). Método de resolución de sistemas de ecuacioneslineales basadas en la utilización de un parámetro de sobrerrelajación para acelerarla convergencia del método iterativo.Q QUICK (Quadratic Upstream Interpolation for Convective Kinetics). Algoritmos deinterpolación cuadrática usados en la discretización de las ecuaciones en losmétodos de volúmenes finitos.R RANS (Reynolds Averaged Navier-Stokes Equations). Ecuaciones de Navier-Stokespromediadas según esfuerzos de Reynolds. Técnicas de solución de las ecuacionesde gobierno discutidas en el capítulo cuatro de este documento. RHS (Right Hand Side). Términos al lado derecho de una igualdad matemática. RIT (Rotor Inlet Temperature). Temperatura de entrada en el rotor. No se trata de untérmino exclusivo de las técnicas numéricas, sino que es un término general usadoen turbomáquinas. 133. Técnicas Numéricas en Mecánica de Fluidos133 RNG (Renormalization Group Theory). Teoría del grupo renormalizado. Técnicamatemática utilizada como modelo de turbulencia. RSM (Reynolds Stresses Model). Modelo de las tensiones de Reynolds. También llamadomodelo de cierre de segundo orden para las ecuaciones de conservación. RSE (Reynolds Stresses Equations). Ecuaciones que utiliza el método definidoanteriormente (RSM). RTE (Radiactive Transfer Equation). Ecuación de transferencia de calor por reacciónradioactiva.S SCBR (Steady Coupled Blade Row). Resolución de las ecuaciones de la Mecánica de Fluidos en turbomáquinas considerando la interacción entre varias coronas de álabes. Modelos estacionarios. SGS (Subgrid Scale). Escala de la malla inferior en un sistema de resolución por medio de varios mallados (multimallado). SIMPLE (Semi Implicit Method for Pressure Linked Equations). Método de resolución de las ecuaciones de gobierno establecido por Patankar (1976). SIMPLEC (Semi Implicit Method for Pressure Linked Equations Consistent). Una variación del método anteriormente descrito. SIP (Strongly Implicit Method). Métodos de resolución de sistemas de ecuaciones lineales basados en algoritmos altamente implícitos. SLIP (Symmetric Limited Possitive Schemes). Algoritmos definidos que conducen a sistemas de ecuaciones simétricos y con coeficientes positivos. SPD (Symmetric Positive Definite Matrix). Matriz simétrica y definida positiva. SSBR (Steady Single Blade Row). Resolución de las ecuaciones en turbomáquinas considerando únicamente el canal entre álabes de una única corona. SUDS (Skew Upwind Differencing Scheme). Método de resolución de las ecuaciones basado en la diferenciación aguas arriba. Establecido por Raithby en 1976. SU/PG (Streamline Upwind Petrov Galerkin). Método espectral basado en la variación de la función peso en elementos finitos para considerar la dirección del flujo. Desarrollado por Brook y Huges en 1982.T TDMA (Tri-Diagonal Matrix Algorithm). Algoritmo de resolución basado en la utilización de sistemas lineales tridiagonales. TVD (Total Variation Diminishing). Algoritmo de resolución utilizado en modelos de flujos compresibles que controla la variación de la solución de forma no lineal. Introducido por Harten en 1983. TSL (Modelo Thin shear Layer). Uno de los modelos de aproximación de las ecuaciones de Navier-Stokes. Se explica más en profundidad en el texto (capítulo cuatro). TKE (Turbulent Kinetic Energy). Energía cinética turbulenta, normalmente nombrada con la letra k. TSTE (Taylor Series Truncation Error). Error de truncatura de los distintos modelos que se basan en los desarrollos en series de Taylor.U UDS (User Defined Scalar). Escalar definido por el usuario. Término usado en los códigos numéricos comerciales. 134. 134 Técnicas Numéricas en Mecánica de Fluidos U-P (Velocity – Pressure). Modelos que usan en su desarrollo las variables primitivas, esdecir, la presión y la velocidad. USLIP (Upwind Symmetric Limited Possitive Schemes). Esquemas aguas arriba ysimétricos con valores positivos.V VOF (Volume Of Fluid). Algoritmo que intenta capturar el movimiento de las partículas siguiendo la fracción de flujo que atraviesa cada cara del mallado. VTE (Vorticity Transport Equations). Ecuaciones de transporte de la turbulencia. Suele darse este nombre a las técnicas numéricas que utilizan la vorticidad.W WSGGM (Weighted Sum of Gray Gases Model). Modelo de la suma de los gases grises. Es un modelo usado en combustión.AIII.2.- Instituciones y organismos que aparecen habitualmente en textosnuméricos. AGARD (Advisory Group for Aerospace Research & Development). Consejo para la investigación y desarrollo aeroespacial. Perteneciente a la OTAN. AIAA (American Institute of Aeronautics and Astronautics). Instituto Americano par la aeronáutica y la astronáutica. API (American Petroleum Institute). Instituto Americano del petróleo. Estados Unidos. ASME (American Society of Mechanical Engineers). Asociación Americana de Ingenieros Mecánicos. Entre otras muchas actividades, publican los PTC (Performanes Test Code) o códigos para ensayos experimentales. ASPE (American Society of Plumbing Engineers). Sociedad Americana de Ingenieros dedicados al diseño e instalación de conductos. Estados Unidos. ASRHAE (American Society of Heating, Refrigeration, and Air Conditioning Engineers). Sociedad Americana de Ingenieros dedicados a instalaciones de calefacción, refrigeración y aire acondicionado. BHRA (British Hidromechanics Research Association). Asociación Británica para la investigación hidromecánica. Inglaterra. BRITE (Basic Research In Technology for Europe). Programa Europeo para el desarrollo de investigación tecnológica básica. CALTECH (California Institute of Technology). Instituto de Tecnologías en California (Estados Unidos). CETIM (Centre Technique d’Industries Mécaniques). Centro Tecnológico de Industrias Mecánicas. Situado en Nantes, Francia. CGPM (Conférence Général des Poids et Mesures). Comité Internacional de Pesas y medidas. Creado en Francia, estableció en 1960 el SI (Sistema Internacional) de medidas. CREMHYG (Centre de Recherche et d’Essais de Machines Hydrauliques de Grenoble). Centro de Investigación y ensayo de máquinas hidráulicas de Grenoble, Francia. DFVLR (Deutsche Forschungs-und Versuchsanstalt für Luft- und Raumfahrt). Laboratorio para la investigación en Alemania. EUROMECH (European Mechanics Society). Sociedad Europea para la Ingeniería Mecánica. 135. Técnicas Numéricas en Mecánica de Fluidos135 GALCIT (Graduate Aeronautical Laboratories at California Institute of Technology). Laboratorios de aeronáutica en el CALTECH. HENSA (Higher Education National Software Archive). Archivo Nacional de Software dedicado a la Educación. Institución de Inglaterra. HVAC (Automotive Heating, Ventilation and Air Conditioning engineers). Asociación de los Ingenieros dedicados a la calefacción, ventilación y aire acondicionado en automóviles. Estados Unidos. IGTI (International Gas Turbine Institute). Instituto Internacional dedicado a las Turbinas de Gas. Organismo dependiente del ASME. ISABE (International Society for Air Breathing Engines). Sociedad Internacional para los Motores de Combustión que utilizan el aire. JPL (Jet Propulsion Laboratory). Es una de la instituciones o laboratorios que funcionan dentro del CALTECH. JSME (Japanese Society of Mechanical Engineers). Sociedad Japonesa de Ingenieros Mecánicos. LLNL (Lawrence Livermore National Laboratory). Laboratorio Americano dedicado a la Mecánica de Fluidos. Estados Unidos. MIT (Massachussets Institute of Technology). Instituto Tecnológico de Massachussets. Uno de los grandes centros de Investigación a nivel mundial. MTU (Motoroen-und Turbinen-Union Munchen). Instituto para la investigación en Motores y Turbinas. Munich (Alemania). NACA (National Advisory Committee for Aeronautics). Comité Nacional para la Aeronáutica. Es el organismo predecesor de la NASA. Estados Unidos. NASA (National Aeronautics and Space Administration). Administración Nacional para la Aeronáutica y el Espacio. Estados Unidos. NGTE (National Gas Turbine Establisment). Instituto para las Turbinas de Gas. Estados Unidos. NIST National Institute of Standards and Technology). Instituto Nacional para la Normativa y las Tecnologías. Estados Unidos. NRC (Nuclear Regulatory Commission). Comisión Reguladora de actividades Nucleares. Estados Unidos. ONERA (Office National d’Etudes et de Recherche Aéroespatiales). Oficina Nacional para la Investigación. Paris, Francia. RAE (Royal Aircraft Establishment). Academia Real para la Aviación. Inglaterra. ROMAC (Rotating Machinery and Controls). Máquinaria rotativa y controles. Centro de investigación de la Universidad de Virginia (Estados Unidos). SHF (Société Hydrothermique de France). Sociedad para estudios Hidráulicos y Térmicos de Francia. UMIST (University of Manchester Institute of Science and Technology). Instituto de Ciencia y Tecnología de la Universidad de Manchester (Inglaterra). UTRC (United Technologies Research Center). Centro de investigación en varias áreas de la ingeniería. Connecticut (Estados Unidos). VKI (von Karman Institute). Instituto de Investigación von Karman. Dedicado a investigación sobre Mecánica de Fluidos. Situado en Bélgica. 136. 136 Técnicas Numéricas en Mecánica de Fluidos 137. ISBN 84-607-9546-2