Aplicaciones del método de MacCormack a diversos problemas fluidomecánicos

88 Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Aplicaciones del método de MacCormack a diversos problemas fluidomecánicos Autor: Manuel Carreño Ruiz Tutor: Miguel Pérez-Saborid Sánchez-Pastor Dep. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2016 2 3 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Aplicaciones del método de MacCormack a diversos problemas fluidomecánicos Autor: Manuel Carreño Ruiz Tutor: Miguel Pérez-Saborid Sánchez-Pastor Profesor titular Dep.Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2016 4 5 Trabajo Fin de rado: Aplicaciones del método de MacCormack a diversos problemas fluidomecánicos Autor: Manuel Carreño Ruiz Tutor: Miguel Pérez-Saborid Sánchez-Pastor El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2016 6 El Secretario del Tribunal 7 A mi familia, amigos y compañeros. 8 Agradecimientos Deseo expresar mi más profundo agradecimiento a todos los docentes con los que he tenido el placer de recorrer este duro camino que ahora cierro con este documento, especialmente a Don Miguel Pérez-Saborid Sánchez-Pastor por el inestimable apoyo recibido durante estos meses de trabajo. También me gustaría agradecer a mi familia y amigos el apoyo que me han brindado durante estos cuatro años, sin ningún asomo de duda sin ellos hoy no estaría aquí. Por último, me gustaría agradecer a mis compañeros la ayuda recibida, especialmente a María Paz Luque Alcántara y Carmen Velarde López de Ayala, cuya compañía durante largos días de trabajo han hecho que todo sea mucho más sencillo y han contribuido de manera inigualable a mi formación tanto personal como universitaria. Manuel Carreño Ruiz Sevilla, 2016 9 Resumen Estas memorias recogen unas instrucciones y pasos sencillos mediante los cuales el lector pueda ser capaz de tener una comprensión del método de MacCormack que le permitirá implementarlo en un software de cálculo numérico como MATLAB o similar para resolver diversos problemas fluidodinámicos. Además de anlizar la validez del método comparándo con resultados de la literatura se procederá a analizar los resultados de problemas como el flujo alrededor de un cilindro cuadrado, la convección en una cavidad cuadrada llevada por la tapadera o el problema de convección natural de Rayleigh-Bénard. Los dos problemas estudiados en más profundidad son la emisión de torbellinos que se produce a partir un número de Reynolds cuando una corriente rebordea un obstáculo, en este caso un cuadrado, y algunos aspectos relacionados con el problema de convección de Rayleigh-Benard, como el número de Reynolds crítico o la influencia de las condiciones de contorno en las paredes de la cavidad. Como aspecto diferenciador a la mayoría de resultados de la bibliografía para el problema de RayleighBénard en esta memoria se contempla la resolución de las ecuaciones de Navier-Stokes compresibles, sin aplicar la aproximación de Boussinesq lo que nos permirá comprobar la validez de esta. 16 Henri Bénard llevó a cabo una serie de experimentos en capas de fluidos de elevada conductividad térmica calentadas desde abajo y cuya su superficie superior estaba expuesta al aire del ambiente, lo que dio comienzo al estudio de los fenómenos de convección natural. Bérnard pudo observar la existencia de un movimiento en el fluido que alcanzaba un régimen estacionario, tras un régimen transitorio inicial, donde se distinguía en la superficie una estructura similiar a la de un panal de abejas de celdas hexagonales regulares como se observa en la figura 1-4. Realizó diversas investigaciones del flujo en el interior de la celda, las cuales revelaron que el fluido ascendía por el centro de la celda y descendía a lo largo del perímetro hexagonal. Figura( 1-4 ) El primer investigador que desarrolló una teoría que ponía de manifiesto de forma clara los mecanismos físicos involucrados en el fenómeno de la convección natural, fue Lord Rayleigh. Rayleigh, resaltó que además de las fuerzas de flotabilidad, deben tenerse en cuenta, cuando un fluido cuya densidad no es uniforme se encuentra en presencia de la gravedad; los mecanismos que tienden a contrarrestar su efecto: la fricción, originada por las fuerzas de viscosidad, y el de conducción de calor, que tiende a homogeneizar en el seno del fluido el campo de temperaturas y, por tanto, también el de densidades. Rayleigh consideró una capa de fluido entre dos placas paralelas infinitas a diferente temperatura dispuestas perpendicularmente a la dirección de la gravedad, con la placa inferior a mayor temperatrura. Rayleigh linealizo las ecuaciones de Navier-Stokes entorno a la situación de equilibrio, descubriendo así que en el sistema lineal de ecuaciones en derivadas parciales aparece el parámetro adimensional, denominado número de Rayleigh: 𝑅𝑎=𝑔𝛽Δ𝑇𝐻3 𝜈𝛼 ( 1-1) donde β es el coeficiente de dilatación térmica, ν la viscosidad cinemática, α la difusividad térmica del fluido, g es la aceleración de la gravitatoria, H es una longitud característica del problema y ΔT es la diferencia de temperatura entre ambas placas. El número de Rayleigh es el Número adimensional que mide la importancia relativa entre los efectos de las fuerzas de flotabilidad y los efectos viscosos y de conducción térmica. Existe un cierto valor umbral del número de Rayleigh a partir del cual los efectos de flotabilidad dominan sobre los de conducción de calor y los viscosos y se desarrolla la convección natural, este se denomina número de Rayleigh crítico. El valor teórico de 𝑅𝑎𝑐=1708, que coincide muy bien con los valores medidos, alrededor de 1708±50, para el caso de dos placas infinitas. Muchos investigadores realizaron mediciones de este número de Rayleigh Crítico para cavidades finitas con diferentes relaciones de aspecto y distintas condiciones de contorno en las paredes, por tanto, existe una gran variedad de bibliografía para comparar los resultados. Finalmente comentar que recientemente se descubrió que los patrones convectivos observados por Bénard no eran debidos a la diferencia de tensión superficial que genera la diferencia de temperaturas, no obstante, este fenómeno sigue denominándose comúnmente como convección de Rayleigh-Bénard. 17 2 ESTRUCTURA DEL TRABAJO En el presente trabajo se muestra un estudio del método numérico explícito de MacCormack desarrollado por R.W.MacCormack y su aplicación a diversos problemas fluidodinámicos para hacer posible su resolución mediante un programa MATLAB relativamente sencillo. En el capítulo 3 analizaremos los fundamentos de dicho método, incluyendo estrategias que hacen que el método sea más estable. Se presentarán las ecuaciones de Navier-Stokes en forma conservativa, esta disposición de las ecuaciones es esencial para poder aplicar este esquema numérico de forma eficiente. Además se incluye una pequeña aplicación de un caso unidimensional para ayudar a entender como funciona la mecánica del método antes de comenzar a resolver problemas más complejos y un ejemplo de un problema bidimensional para comenzar a entender el método aplicado a las ecuaciones de Navier-Stokes. Un caso muy simple, una caja cuadrada llena de fluido cuya tapadera se desplaza con velocidad U constante. Se simplifican las ecuaciones de Navier-Stokes completas mostradas al comienzo del apartado y se aplican las condiciones de contorno correspondientes. Este ejemplo nos sirve como verificador del método ya que ha sido resuelto por muchos autores y nos permite comparar los resultados. El capítulo 4, es una extensión del apartado 4 del capítulo 3, una cavidad rectangular alargada, que se puede asimilar a un conducto, con un cuadrado en medio que hace sus veces de obstáculo para el fluido. Se imponen las condiciones de contorno correspondientes tal que el problema sea equivalente al cuadrado viajando a velocidad constante en el conducto. El objetivo de esta aplicación será el cálculo de la resistencia generada por el cuadrado asi como la sustentación si la hubiese. También vamos a estudiar la sensibilidad a la relación de aspecto de la cavidad. Una vez llevado a cabo este análisis podremos elegir una relación de aspecto de la cavidad de tal forma que el resultado se aproxime a el cuadrado inmerso en una corriente uniforme y entonces analizaremos el problema de ‘Vortex shedding’ que comentamos anteriormente, calculando la frecuencia con la que se emiten torbellinos aguas abajo. En el capítulo 5 se contempla la resolución de un problema parecido al anterior pero en este caso vamos a incluir el problema térmico. El cuadrado estará a una temperatura superior a la del fluido por tanto el fluido actuará como refrigerante llevándose el calor que pudiese desprender el cuadrado. Se analizarán los perfiles de temperatura a distintos números de Prandtl. El capítulo 6, recoge un estudio del fenómeno de inestabilidad térmica que se produce cuando una cavidad es calentada por su parte inferior en presencia de las fuerzas másicas. Esta inestabilidad produce la conocida como convección natural o convección de Rayleigh-Bénard. Este ejemplo es algo más complejo y por tanto se desarrollará un breve tratamiento teórico del problema. Se realizarán estudios en relación a las distintas celdas convectiva que se forman al variar la relación de aspecto y la influencia de los números adimensionales de nuestro problema en la aparición y vigorosidad de la convección obtenida. Estos análisis se han realizado sin hacer uso de la aproximación de Boussinesq, que es lo usual en la literatura. Esto permitirá ver el efecto de la compresibilidad en el régimen de convección natural. 18 3 MÉTODO DE MACCORMACK 3.1 Elección del método Existen infinidad de métodos numéricos para resolver problemas en mecánica de fluidos, cada uno tiene sus ventajas y sus inconvenientes. Según R.W. MacCormack [28] todo principiante en CFD comienza utilizando métodos de diferencias finitas por la familiaridad con las ecuaciones diferenciales y después paulatinamente tienden a métodos basados en volúmenes finitos. Como el objetivo primordial de este documento es la aplicación de un método numérico a problemas de convección para que pueda ser utilizado por alumnos de asignaturas de mecánica de fluidos hemos creido mas apropiado comenzar con un método de diferencias finitas. Existen muchos métodos de diferencias finitas como los métodos de Lax-Friedrichs, Lax-Wendroff o MacCormack entre otros. Hemos escogido el método de MacCormack por que además de proporcionar resultados más que satisfactorios para diversos problemas fluidodinámicos, es el más sencillo de programar y por tanto una muy buena elección para iniciarse en la mecánica de fluidos computacional. También el Método de MacCormack admite una extensión a problemas tridimensionales bastante sencilla lo cual deja abierta una continuación de este trabajo para resolver problemas fluidodinámicos tridimensionales. 3.2 FundamentoTeórico Este método consiste en dos etapas, una de predicción y otra correctora. La primera consiste en una aproximación de las derivadas espaciales mediante la diferencia progresiva y la segunda mediante la diferencia regresiva. Después se realiza una media aritmética entre ambas generando una aproximación de segundo orden equivalente a la diferencia central como se demostrará en el siguiente apartado. Las derivadas temporales se aproximan mediante la diferencia progresiva. Véase el siguiente esquema general: 𝐸𝑐𝑢𝑎𝑐𝑖ó𝑛: 𝜕(𝑈) 𝜕𝑡 +𝜕(𝐹) 𝜕𝑥 +𝜕(𝐺) 𝜕𝑦 +𝑆=0 ( 3-1) 𝑃𝑟𝑒𝑑𝑖𝑐𝑡𝑜𝑟:𝑈𝑖,𝑗 𝑛+1 =𝑈𝑖,𝑗 𝑛−𝛥𝑡(𝐹𝑖+1,𝑗 𝑛−𝐹𝑖,𝑗 𝑛 𝛥𝑥 )−𝛥𝑡(𝐺𝑖,𝑗+1 𝑛−𝐺𝑖,𝑗 𝑛 𝛥𝑦 )−𝛥𝑡𝑆 ( 3-2) 𝐶𝑜𝑟𝑟𝑒𝑐𝑡𝑜𝑟:𝑈𝑖,𝑗 𝑛+1=12(𝑈𝑖,𝑗 𝑛+1 +𝑈𝑖,𝑗 𝑛−𝛥𝑡(𝐹𝑖,𝑗 𝑛+1 −𝐹𝑖−1,𝑗 𝑛+1  𝛥𝑥 )−𝛥𝑡(𝐺𝑖,𝑗 𝑛+1 −𝐺𝑖,𝑗−1 𝑛+1  𝛥𝑦 )−𝛥𝑡𝑆) ( 3-3) Como se muestra en el esquema anterior la ecuación debe de transformarse en forma conservativa para la ejecución de este método. Se puede intercambiar el orden de las aproximaciones regresiva y progresiva en el esquema superior, es más, es recomendable intercambiarlos alternativamente para evitar direcciones espaciales preferentes. En la siguiente tabla se muestran las posibles combinaciones: 19 Tabla ( 3-1 ) Predictor Corrector Dirección-x Dirección-y Dirección-x Dirección-y DP DP DR DR DB DR DP DP DP DR DR DP DR DP DP DR En muchos casos alternar estos ciclos de derivadas proporciona una mejora en la estabilidad numérica. No existe una demostración para la convergencia del método en problemas bidimensionales, pero la siguiente condición se ha mostrado efectiva en la práctica como veremos más adelante: 𝛥𝑡≤1 |𝑢| 𝛥𝑥+|𝑣| 𝛥𝑦+𝑐√1 𝛥𝑥2+1 𝛥𝑦2 ( 3-4) Se conoce como condición de Courant-Friedrich-Levy (CFL) Es interesante notar que este método puede originar un falso estado Pseudo-transitorio al alcanzar la unidad de redondeo de la computadora y no terminar de converger debido a sus dos etapas. Como se puede apreciar en la ecuación anterior, para la aplicación de este método es esencial que las ecuaciones estén en forma conservativa. A continuación se muestran las ecuaciones de Navier-Stokes completas para el problema bidimensional en esta forma: 𝜕(𝑈) 𝜕𝑡 +𝜕(𝐹) 𝜕𝑥 +𝜕(𝐺) 𝜕𝑦 +𝑆=0 ( 3-5) 𝑈=[𝜌 𝜌𝑢 𝜌𝑣 𝐸]=[𝑈1 𝑈2 𝑈3 𝑈4] ( 3-6) 𝐹= [ 𝜌𝑢 𝜌𝑢2+𝑝−𝜏𝑥𝑥 𝜌𝑢𝑣−𝜏𝑥𝑦 (𝐸+𝑝)𝑢−𝑢𝜏𝑥𝑥−𝑣𝜏𝑥𝑦+𝑞𝑥 ] ( 3-7) 20 𝐺= [ 𝜌𝑣 𝜌𝑢𝑣−𝜏𝑥𝑦 𝜌𝑣2+𝑝−𝜏𝑦𝑦 (𝐸+𝑝)𝑣−𝑢𝜏𝑥𝑦−𝑣𝜏𝑦𝑦+𝑞𝑦 ] ( 3-8) El vector S, que no contiene derivadas espaciales ni temporales, se denomina término fuente y se reserva para las fuerzas externas que puedan aparecer en el fluido como pueden ser la gravedad o adición de calor. Por ejemplo: 𝑆=[00 𝜌𝑔 𝜌𝑣𝑔] ( 3-9) Una vez calculados los valores de U se pueden ‘decodificar’ las variables primarias de la siguiente forma: 𝜌=𝑈1 𝑢=𝑈2 𝑈1 𝑣=𝑈3 𝑈1 ( 3-10) 𝑒=𝑈4 𝑈1−(𝑈2 𝑈1)2+(𝑈3 𝑈1)2 2 ( 3-11) Que junto a las ecuaciones de estado nos permitirán determinar todas las variables fluido-dinámicas. En este manual se considerará el caso de un gas perfecto: 𝑒=𝑐𝑣𝑇 𝑝=𝜌𝑅𝑔𝑇 ( 3-12) 3.3 Aproximación numérica en mallas equiespacidas El esquema de MacCormack consiste en utilizar una etapa de diferencia progresiva y otra etapa de diferencia regresiva y realizar una combinación entre ambas. Estas estimaciones son cada una de primer orden pero la combinación de ambas es equivalente a una diferencia central que cómo veremos a continuación es de Segundo orden. Las aproximaciones más comunes para derivadas de primer orden son: 21  Diferencia progresiva: 𝑑𝜙 𝑑𝑥≈𝜙𝑖+1−𝜙𝑖 Δ𝑥 ( 3-13)  Diferencia regresiva: 𝑑𝜙 𝑑𝑥≈𝜙𝑖−𝜙𝑖−1 Δ𝑥 ( 3-14)  Diferencia central: 𝑑𝜙 𝑑𝑥≈𝜙𝑖+1−𝜙𝑖−1 2Δ𝑥 ( 3-15) Si se realiza una expansión en serie de Taylor de la función 𝜙 alrededor del punto 𝑥𝑖 obtendremos: 𝜙𝑖+1=𝜙𝑖+𝛥𝑥𝑑𝜙 𝑑𝑥𝑖+𝛥𝑥2 2𝑑2𝜙 𝑑𝑥2𝑖+𝛥𝑥3 6𝑑3𝜙 𝑑𝑥3𝑖+⋯ ( 3-16) 𝜙𝑖=𝜙𝑖 ( 3-17) 𝜙𝑖+1=𝜙𝑖−𝛥𝑥𝑑𝜙 𝑑𝑥𝑖+𝛥𝑥2 2𝑑2𝜙 𝑑𝑥2𝑖−𝛥𝑥3 6𝑑3𝜙 𝑑𝑥3𝑖+⋯ ( 3-18) Substituyendo estos desarrollos en las aproximaciones anteriores tenemos que:  Diferencia progresiva: 𝜙𝑖+1−𝜙𝑖 Δ𝑥 =𝑑𝜙 𝑑𝑥𝑖+Δ𝑥 2𝑑2𝜙 𝑑𝑥2𝑖+⋯ ( 3-19) El error de truncamientoes de primer orden.  Diferencia regresiva: 𝑑𝜙 𝑑𝑥=𝜙𝑖−𝜙𝑖−1 Δ𝑥 ( 3-20) 𝜙𝑖−𝜙𝑖−1 Δ𝑥 =𝑑𝜙 𝑑𝑥𝑖−Δ𝑥 2𝑑2𝜙 𝑑𝑥2𝑖+⋯ ( 3-21) El error de truncamiento es de primer orden. 22  Diferencia central: 𝜙𝑖+1−𝜙𝑖−1 2𝛥𝑥 =𝑑𝜙 𝑑𝑥𝑖−𝛥𝑥2 6𝑑3𝜙 𝑑𝑥3𝑖+⋯ ( 3-22) El error de truncamientoes de Segundo orden. Es interesante comentar que si sumamos la diferencia progresiva y la diferencia regresiva y dividios entre 2, que es esencialmente lo que hace el método de MacCormack, obtenemos la diferencia central que como se acaba de comprobar es de segundo orden. Para actualizar las condiciones de contorno en algunos casos necesitaremos aplicar diferencias progresivas y regresivas de segundo orden:  Diferencia progresiva: 𝑑𝜙 𝑑𝑥𝑖=−3𝜙𝑖−4𝜙𝑖+1+𝜙𝑖+2 2Δ𝑥 +𝑂(Δ𝑥2) ( 3-23)  Diferencia regresiva: 𝑑𝜙 𝑑𝑥𝑖=3𝜙𝑖−4𝜙𝑖−1+𝜙𝑖−2 2Δ𝑥 +𝑂(Δ𝑥2) ( 3-24) En algunos casos es necesario realizar aproximaciones de derivdas de orden superior, por ejemplo, para actualizar condiciones de contorno o para computar los términos correspondientes a efectos de viscosidad. Las más relevantes son:  Diferencia central: 𝜙𝑖+1−2𝜙𝑖+𝜙𝑖−1 Δ𝑥2 ( 3-25)  Diferencia progresiva: 𝜙𝑖+2−2𝜙𝑖+1+𝜙𝑖 Δ𝑥2 ( 3-26)  Diferencia regresiva: 𝜙𝑖−2𝜙𝑖−1+𝜙𝑖−2 Δ𝑥2 ( 3-27) Que con un procedmientocompetamente análogo al realizado en el apartado anterior se puede demostrar que la diferencia central es de segundo orden de precisión y las otras dos diferencias únicamente de primer orden. En algunas ocasiones será preciso aproximar derivadas de segundo orden con precisión de orden dos lo que realizaremos de la siguiente forma: 23  Diferencia progresiva: 𝑑2𝜙 𝑑𝑥2𝑖=−2𝜙𝑖−5𝜙𝑖+1+4𝑢𝑖+2−𝜙𝑖+3 Δ𝑥2 ( 3-28)  Diferencia regresiva: 𝑑2𝜙 𝑑𝑥2𝑖=2𝜙𝑖−5𝜙𝑖−1+4𝑢𝑖−2−𝜙𝑖−3 Δ𝑥2 ( 3-29)  Diferencia cruzada progresiva: 𝜕2𝜙 𝜕𝑥𝜕𝑦𝑖,𝑗=−−(𝜙𝑖+2,𝑗−1−𝜙𝑖+2,𝑗+1)+4(𝜙𝑖+1,𝑗−1−𝜙𝑖+1,𝑗+1)−3(𝜙𝑖,𝑗−1−𝜙𝑖,𝑗+1) 4ΔxΔy +𝑂(Δ𝑥2,Δ𝑦2,Δ𝑥Δ𝑦) ( 3-30)  Diferencia cruzada regresiva: 𝜕2𝜙 𝜕𝑥𝜕𝑦𝑖,𝑗=−−(𝜙𝑖−2,𝑗+1−𝜙𝑖−2,𝑗−1)+4(𝜙𝑖−1,𝑗+1−𝜙𝑖−1,𝑗−1)−3(𝜙𝑖,𝑗+1−𝜙𝑖,𝑗−1) 4ΔxΔy +𝑂(Δ𝑥2,Δ𝑦2,Δ𝑥Δ𝑦) ( 3-31) Con estas sencillas fórmulas podremos aplicar las condiciones de contorno y las aproximaciones que requieren las ecuaciones. Las ecuaciones diferenciales que se muestran en el apartado anterior han de ser discretizadas para poder aplicar este método de diferencias finitas. Se ha escogido una malla de (𝑁𝑥+1)×(𝑁𝑦+1) nodos donde se evaluarán las magnitudes fluidodinámicas necesarias. La distancia entre cada nodo se denominará Δ𝑥 en dirección x y Δ𝑦 en dirección y. Vamos a diferenciar entre nodos interiores y nodos frontera, los nodos interiores (2:𝑁𝑥)×(2:𝑁𝑦) que será donde se aplicarán y calcularan las magnitudes fluidodinámicas provenientes de las ecuaciones de NavierStokes y los nodos frontera donde será necesario imponer condiciones de contorno, como veremos más adelante. Es fundamental aplicar las ecuaciones en forma compacta ya que la condición CFL nos limita mucho el intervalo temporal máximo que podemos escoger generando unos tiempos de computación considerablemente elevados. Por tanto la incorporación de bucles nos demoraría aún mas lo cual produciría tiempos de computación inabordables. A continuación mostramos como se han de implementar las ecuaciones, utilizaremos la ecuación de continuidad a modo de ejemplo. 𝜕𝜌 𝜕𝑡=−𝜕(𝜌𝑢) 𝜕𝑥 −𝜕𝜌𝑣 𝜕𝑦 ( 3-32) 24 Un aspecto clave que simplifica mucho la implementación es la definición de variables auxiliares sobre las cuales aplicamos directamente la derivada, por ejemplo, en la ecuación de continuidad estas variables serían 𝜌𝑢 𝑦 𝜌𝑣 que pueden ser calculadas trivialmente a partir de la velocidad y la densidad. En primer lugar vamos a considerar el término predictor, en esta etapa utilizábamos diferencia progresiva. Recordemos que tanto la densidad como el resto de las variables fluidodinámicas dependen del tiempo y de la posición en el espacio, por tanto los términos de la ecuación quedarían como sigue. 𝜕𝜌 𝜕𝑡≅𝜌𝑖,𝑗 𝑡+1−𝜌𝑖,𝑗𝑡 Δ𝑡 ( 3-33) 𝜕𝜌𝑢 𝜕𝑥≅𝜌𝑢𝑖+1,𝑗 𝑡−𝜌𝑢𝑖,𝑗𝑡 Δ𝑥 ( 3-34) 𝜕𝜌𝑣 𝜕𝑦≅𝜌𝑣𝑖,𝑗+1 𝑡−𝜌𝑣𝑖,𝑗𝑡 Δ𝑦 ( 3-35) Quedando la ecuación: 𝜌𝑖,𝑗 𝑡+1−𝜌𝑖,𝑗𝑡 Δ𝑡 =−𝜌𝑢𝑖+1,𝑗 𝑡−𝜌𝑢𝑖,𝑗𝑡 Δ𝑥 −𝜌𝑣𝑖,𝑗+1 𝑡−𝜌𝑣𝑖,𝑗𝑡 Δ𝑦 ( 3-36) En la cual observamos que todos los valores son conocidos en el instante t siendo la única incognita la densidad en el instante t+1 que podemos calcular como: 𝜌𝑖,𝑗 𝑡+1=𝜌𝑖,𝑗𝑡−Δ𝑡𝜌𝑢𝑖+1,𝑗 𝑡−𝜌𝑢𝑖,𝑗𝑡 Δ𝑥 −Δ𝑡𝜌𝑣𝑖,𝑗+1 𝑡−𝜌𝑣𝑖,𝑗𝑡 Δ𝑦 ( 3-37) Que en forma matricial se puede expresar como: 𝜌(2:𝑁𝑥,2:𝑁𝑦)𝑡+1 =𝜌(2:𝑁𝑥,2:𝑁𝑦)𝑡−Δ𝑡 Δ𝑥(𝜌𝑢(3:𝑁𝑥+1,2:𝑁𝑥)𝑡−𝜌𝑢(2:𝑁𝑥,2:𝑁𝑥)𝑡) −Δ𝑡 Δ𝑦(𝜌𝑣(2:𝑁𝑥,3:𝑁𝑥+1)𝑡−𝜌𝑣(2:𝑁𝑥,2:𝑁𝑥)𝑡) ( 3-38) A este valor lo renombraremos como 𝜌𝑖,𝑗 𝑡+1       por ser el término predictor. Ahora realizaremos la misma operación pero para el término corrector, en el que recordemos que se aproximaban sus términos mediante la diferencia regresiva. 𝜕𝜌 𝜕𝑡≅𝜌𝑖,𝑗 𝑡+1−𝜌𝑖,𝑗𝑡+1       Δ𝑡 ( 3-39) 25 𝜕𝜌𝑢 𝜕𝑥≅𝜌𝑢𝑖,𝑗 𝑡+1       −𝜌𝑢𝑖−1,𝑗𝑡+1       Δ𝑥 ( 3-40) 𝜕𝜌𝑣 𝜕𝑦≅𝜌𝑣𝑖,𝑗 𝑡+1       −𝜌𝑢𝑖−1,𝑗𝑡+1       Δ𝑦 ( 3-41) Quedando la ecuación: 𝜌𝑖,𝑗 𝑡+1−𝜌𝑖,𝑗𝑡+1       Δ𝑡 =−𝜌𝑢𝑖,𝑗 𝑡+1       −𝜌𝑢𝑖−1,𝑗𝑡+1       Δ𝑥 −𝜌𝑣𝑖,𝑗 𝑡+1       −𝜌𝑢𝑖−1,𝑗𝑡+1       Δ𝑦 ( 3-42) En la cual observamos que todos los valores son conocidos en el instante 𝑡+1        siendo la única incognita la densidad en el instante t+1 que podemos calcular como: 𝜌𝑖,𝑗 𝑡+1=𝜌𝑖,𝑗𝑡+1       −Δ𝑡𝜌𝑢𝑖,𝑗 𝑡+1       −𝜌𝑢𝑖−1,𝑗𝑡+1       Δ𝑥 −Δ𝑡𝜌𝑣𝑖,𝑗 𝑡+1       −𝜌𝑢𝑖−1,𝑗𝑡+1       Δ𝑦 ( 3-43) A este valor de la densidad lo renombraremos como 𝜌𝑖,𝑗 𝑡+1∗ por ser el corrector. 𝜌𝑖,𝑗 𝑡+1=𝜌𝑖,𝑗 𝑡+1∗+𝜌𝑖,𝑗 𝑡 2 ( 3-44) Que en forma matricial se puede expresar como: 𝜌(2:𝑁𝑥,2:𝑁𝑦)𝑡+1 =0.5(𝜌(2:𝑁𝑥,2:𝑁𝑦)𝑡+𝜌(2:𝑁𝑥,2:𝑁𝑦)𝑡+1       −Δ𝑡 Δ𝑥(𝜌𝑢(2:𝑁𝑥,2:𝑁𝑥)𝑡+1       −𝜌𝑢(1:𝑁𝑥−1,2:𝑁𝑥)𝑡+1       ) −Δ𝑡 Δ𝑦(𝜌𝑣(2:𝑁𝑥,2:𝑁𝑥)𝑡+1       −𝜌𝑣(2:𝑁𝑥,1:𝑁𝑥−1)𝑡+1       )) ( 3-45) Aquí únicamente se ha mostrado el desarrollo de la ecuación de continuidad por ser común en todos los ejemplos, en cada aplicación encontrarán un código MATLAB con todas las ecuaciones implementadas. Antes de comenzar problemas más complejos es conveniente ver un ejemplo unidimensional para captar la esencia del método. 32 %====================ECUACIÓN DE ONDA====================% %LAS VARIABLES P_ SE CORRESPONDEN A LOS TÉRMINOS DEL PREDICTOR% %========================================================% close all;clear all; nx=0;%NÚMERO DE INTERVALOS DE LA PARTICIÓN T=5;%%TIEMPO TRANSCURRIDO c=1;%%CONSTANTE DE LA ECUACIÓN DE ONDA Dx=2*pi/nx;%%NORMA DE LA PARTICIÓN ESPACIAL Dt=0.1;%1/c*Dx*0.9;%%NORMA DE LA PARTICIÓN TEMPORAL u(1:nx+1)=sin(0:Dx:2*pi);%CONDICIÓN INICIAL P_u(1:nx+1)=sin(0:Dx:2*pi);%CONDICIÓN INICIAL %==============BUCLE TEMPORAL===========================% for t=1:T/Dt P_u(2:nx)=u(2:nx)-(Dt/Dx)*c*(u(3:nx+1)-u(2:nx));%%PREDICTOR P_u(1)=u(1)-(Dt/2/Dx)*c*(-P_u(3)+4*P_u(2)-3*P_u(1));%%C.C x=0 P_u(nx+1)=u(nx+1)+(Dt/2/Dx)*c*(-u(nx-1)+4*u(nx)-3*u(nx+1));%%C.C x=2PI u(2:nx)=(u(2:nx)+P_u(2:nx)-(Dt/Dx)*c*(P_u(2:nx)-P_u(1:nx1)))*0.5;%%CORRECTOR u(1)=-sin(t*Dt);%%C.C x=0 u(nx+1)=(u(nx+1)+P_u(nx+1)+(Dt/2/Dx)*c*(-P_u(nx-1)+4*P_u(nx)- 3*P_u(nx+1)))*0.5;%%C.C x=2PI end plot(linspace(0,2*pi,length(u)),u,'*')%%PLOT SOLUCIÓN NUMÉRICA hold on p=sin(-t*Dt*c:2*pi/1000:2*pi-t*Dt*c);%%SOLUCIÓN ANALÍTICA plot((linspace(0,2*pi,length(p))),p,'r');%%PLOT SOLUCIÓN ANALÍTICA legend('Solución Numérica','Solución Analítica') title('N_x=50 T=5s') xlabel('x') ylabel('u') %========================================================% 33 3.4 Aplicación: Cavidad cuadrada, con la tapadera desplazándose a velocidad constante. Una vez analizado y comprendido el problema unidimensional vamos a extender el método a dos dimensiones. El problema que se muestra a continuación viene siendo usado como validación de códigos de CFD ya que posee una geometría simple, las condiciones de contorno también lo son y lo más importante, ha sido reproducido por varios autores y existe bibliografía con la que comparar. Figura( 3-7 ) A continuación, se muestran las ecuaciones de Navier-Stokes en las que se han suprimido los términos correspondientes a las fuerzas másicas debido a que estas son despreciables frente a las fuerzas de inercia para los casos analizados: 𝐹𝑟2=𝑈2 𝑔𝐷≫1 ( 3-65) También vamos a considerar que como el fluido se encuentra inicialmente a temperatura constante y todas las paredes se encuentran a esa temperatura, las variaciones de temperatura serán muy pequeñas si las velocidades son pequeñas, como veremos a continuación pudiendo considerar el caso isotermo, realizando estimaciones de órdenes de magnitud en la ecuación de la entalpía: 𝜌𝑐𝑝𝜕𝑇 𝜕𝑡+𝜌𝑐𝑝𝑣 󰇍 󰇍 .∇𝑇=𝜕𝑝 𝜕𝑡+𝑣 󰇍 󰇍 .∇𝑝+∇.(𝜅∇T)+τ′:∇𝑣 ( 3-66) Si se introducen los siguientes números adimensionales: 𝑅𝑒=𝜌𝑈𝐷 𝜇,𝑃𝑟=𝜈𝛼, 𝑀= 𝑢 𝛾𝑅𝑔𝑇0, ( 3-67) 34 Se obtiene:  Incrementos de temperatura debido a variaciones de presión. Δ𝑇𝑝 𝑇0~𝑈2 𝐶𝑝𝑇0~𝑀2 ( 3-68) Donde se ha tenido en cuenta que: Δ𝑝 ~𝑈2 𝐷  Incrementos de temperatura debido a conducción de calor. Δ𝑇𝑞 𝑇0~𝑇1−𝑇2 𝑇0𝑃𝑟𝑅𝑒 ( 3-69)  Incrementos de temperatura debido a disipación viscosa. Δ𝑇𝑑𝑖𝑠 𝑇0~𝑈2 𝐶𝑝𝑇0𝑅𝑒~𝑀2 𝑅𝑒 ( 3-70) En estas tres expresiones apreciamos como si el número de Mach es pequeño, las temperaturas de las paredes son aproximadamente iguales y el número de Reynolds no es pequeño, la aproximación por el caso isotermo es muy adecuada. Por tanto, las ecuaciones de Navier-Stokes bidimensionales se pueden expresar como: 𝐶𝑜𝑛𝑡𝑖𝑛𝑢𝑖𝑑𝑎𝑑: 𝜕𝜌 𝜕𝑡+𝜕(𝜌𝑢) 𝜕𝑥 +𝜕(𝜌𝑣) 𝜕𝑦 =0 ( 3-71) 𝐶𝑎𝑛𝑡𝑖𝑑𝑎𝑑 𝑑𝑒 𝑚𝑜𝑣𝑖𝑚𝑖𝑒𝑛𝑡𝑜 𝑥:𝜕(𝜌𝑢) 𝜕𝑡 +𝜕(𝜌𝑢2+𝑝−4𝜇 3𝑑𝑢 𝑑𝑥−𝜇3𝑑𝑣 𝑑𝑦) 𝜕𝑥 +𝜕(𝜌𝑢𝑣−𝜇𝑑𝑢 𝑑𝑦) 𝜕𝑦 =0 ( 3-72) 𝐶𝑎𝑛𝑡𝑖𝑑𝑎𝑑 𝑑𝑒 𝑚𝑜𝑣𝑖𝑚𝑖𝑒𝑛𝑡𝑜 𝑦:𝜕(𝜌𝑣) 𝜕𝑡 +𝜕(𝜌𝑢𝑣−𝜇𝑑𝑣 𝑑𝑥) 𝜕𝑥 +𝜕(𝜌𝑣2+𝑝−4𝜇 3𝑑𝑣 𝑑𝑦−𝜇3𝑑𝑢 𝑑𝑥) 𝜕𝑦 =0 ( 3-73) ecuación de estado: p=𝑅𝑔𝑇0𝜌 ( 3-74) Siendo c la velocidad del sonido en el medio. 35 Como hemos visto anteriormente a bajas velocidades la temperatura debería ser aproximadamente constante por tanto estas ecuaciones aproximan el límite isotermo. Es conveniente introducir las siguientes variables adimensionales: 𝑥∗=𝑥𝐷,𝑦∗=𝑦𝐷,𝑢∗=𝑢𝑈,𝑣∗=𝑣𝑈,𝑡∗=𝑈 𝐷𝑡,𝜌∗=𝜌 𝜌0,𝑝∗=𝑝 𝜌0𝑈2 ( 3-75) Siendo D un lado de la cavidad, U la velocidad de la tapa y 𝜌0 una densidad de referencia. A partir de ahora se suprimirá el asterisco de las variables adimensionales por simplicidad. Las ecuaciones resultantes son: 𝐶𝑜𝑛𝑡𝑖𝑛𝑢𝑖𝑑𝑎𝑑: 𝜕𝜌 𝜕𝑡+𝜕(𝜌𝑢) 𝜕𝑥 +𝜕(𝜌𝑣) 𝜕𝑦 =0 ( 3-76) 𝐶𝑎𝑛𝑡𝑖𝑑𝑎𝑑 𝑑𝑒 𝑚𝑜𝑣𝑖𝑚𝑖𝑒𝑛𝑡𝑜 𝑥:𝜕(𝜌𝑢) 𝜕𝑡 +𝜕(𝜌𝑢2+𝑝− 4 3𝑅𝑒𝜕𝑢 𝜕𝑥−1 3𝑅𝑒𝜕𝑣 𝑑𝑦) 𝜕𝑥 +𝑑(𝜌𝑢𝑣−1 𝑅𝑒𝜕𝑢 𝜕𝑦) 𝜕𝑦 =0 ( 3-77) 𝐶𝑎𝑛𝑡𝑖𝑑𝑎𝑑 𝑑𝑒 𝑚𝑜𝑣𝑖𝑚𝑖𝑒𝑛𝑡𝑜 𝑦:𝜕(𝜌𝑣) 𝜕𝑡 +𝑑(𝜌𝑢𝑣−1 𝑅𝑒𝜕𝑣 𝜕𝑥) 𝜕𝑥 +𝜕(𝜌𝑣2+𝑝− 4 3𝑅𝑒𝜕𝑣 𝜕𝑦−1 3𝑅𝑒𝜕𝑢 𝜕𝑥) 𝜕𝑦 =0 ( 3-78) 𝑒𝑐𝑢𝑎𝑐𝑖ó𝑛 𝑑𝑒 𝑒𝑠𝑡𝑎𝑑𝑜: 𝑝= 𝜌 𝛾𝑀2 ( 3-79) Sujeto a las condiciones de contorno: 𝑢(0,𝑦)=0 𝑣(0,𝑦)=0 𝑢(1,𝑦)=0 𝑣(1,𝑦)=0 ( 3-80) 𝑢(𝑥,0)=0 𝑣(𝑥,0)=0 𝑢(𝑥,1)=1 𝑣(𝑥,1)=0 ( 3-81) Un factor importante a la hora de implementar satisfactoriamente el método es imponer correctamente las condiciones de contorno en densidad o equivalentemente en presión. Vamos a realizar esto utilizando la ecuación de continuidad: 𝑇𝑎𝑝𝑎 𝑖𝑛𝑓𝑒𝑟𝑖𝑜𝑟: 𝜕𝜌 𝜕𝑡=−𝜕(𝜌𝑣) 𝜕𝑦 ( 3-82) 36 𝑇𝑎𝑝𝑎 𝑠𝑢𝑝𝑒𝑟𝑖𝑜𝑟: 𝜕𝜌 𝜕𝑡=−𝜕𝜌 𝜕𝑥−𝜕(𝜌𝑣) 𝜕𝑦 ( 3-83) 𝑇𝑎𝑝𝑎 𝑖𝑧𝑞𝑢𝑖𝑒𝑟𝑑𝑎: 𝜕𝜌 𝜕𝑡=−𝜕(𝜌𝑢) 𝜕𝑥 ( 3-84) 𝑇𝑎𝑝𝑎 𝑑𝑒𝑟𝑒𝑐ℎ𝑎: 𝜕𝜌 𝜕𝑡=−𝜕(𝜌𝑢) 𝜕𝑥 ( 3-85) Las ecuaciones de Navier Stokes se encuentran discretizadas en el apartado 2, únicamente tendremos que despreciar las fuerzas másicas y al considerar el caso isotermo podremos obviar la ecuación de la energía. Vamos a analizar en más detalle la discretización de las condiciones de contorno en densidad que son algo más complejas. A partir de estas expresiones podemos calcular la densidad en el contorno de la caja mediante las siguientes aproximaciones numéricas que vimos en el segundo apartado. Para la derivada temporal podremos utilizar la diferencia progresiva. 𝜕𝜌 𝜕𝑡≅𝜌𝑖,𝑗 𝑡+1−𝜌𝑖,𝑗𝑡 Δ𝑡 ( 3-86) Sin embargo no se puede utilizar la diferencia progresiva/regresiva indistintamente para una misma etapa ya que al estar en la frontera estas diferencias requerirían nodos fuera del dominio computacional considerado, en otras palabras excederíamos la dimensión de la matriz. Puesto que no se puede seguir el guión general del método, el cual le confiere una precisión de segundo orden en el tiempo y en el espacio, y queremos mantener dicha precisión debemos aplicar aproximaciones de segundo orden. La diferencia central no se puede utilizar ya que excederíamos nuevamente la dimensión de la matriz por tanto elegiremos la diferencia progresiva o regresiva de tal forma que siempre nos mantengamos dentro del dominio computacional. Los valores de la densidad en las esquinas se corresponderán con los valores extremos de las tapas laterales lo cual nos permite aproximar la derivada espacial de la densidad que aparece en la tapa superior mediante la diferencia central. 𝑇𝑎𝑝𝑎 𝑖𝑧𝑞𝑢𝑖𝑒𝑟𝑑𝑎: 𝜕𝜌𝑢 𝜕𝑥𝑖,𝑗≅3𝜌𝑢𝑖,𝑗−4𝜌𝑢𝑖+1,𝑗+𝜌𝑢𝑖+2,𝑗 2Δ𝑥 ( 3-87) 𝑇𝑎𝑝𝑎 𝑑𝑒𝑟𝑒𝑐ℎ𝑎: 𝜕𝜌𝑢 𝜕𝑥𝑖,𝑗≅−3𝜌𝑢𝑖,𝑗−4𝜌𝑢𝑖−1,𝑗+𝜌𝑢𝑖−2,𝑗 2Δ𝑥 ( 3-88) 𝑇𝑎𝑝𝑎 𝑖𝑛𝑓𝑒𝑟𝑖𝑜𝑟: 𝜕𝜌𝑣 𝜕𝑦𝑖,𝑗≅3𝜌𝑣𝑖,𝑗−4𝜌𝑣𝑖,𝑗+1+𝜌𝑣𝑖,𝑗+2 2Δ𝑥 ( 3-89) 37 𝑇𝑎𝑝𝑎 𝑠𝑢𝑝𝑒𝑟𝑖𝑜𝑟: 𝜕𝜌𝑣 𝜕𝑦𝑖,𝑗≅−3𝜌𝑣𝑖,𝑗−4𝜌𝑣𝑖,𝑗−1+𝜌𝑣𝑖,𝑗−2 2Δ𝑥 ( 3-90) 𝜕𝜌 𝜕𝑥𝑖,𝑗 ≅𝜌𝑖+1,𝑗−𝜌𝑖−1,𝑗 2Δ𝑥 ( 3-91) Con estas expresiones no es difícil calcular los valores de la densidad en el instante 𝑡+1        (𝑐𝑜𝑟𝑟𝑒𝑠𝑝𝑜𝑛𝑑𝑖𝑒𝑛𝑡𝑒 𝑎𝑙 𝑝𝑟𝑒𝑑𝑖𝑐𝑡𝑜𝑟) en la frontera del dominio. 𝑇𝑎𝑝𝑎 𝑖𝑧𝑞𝑢𝑖𝑒𝑟𝑑𝑎: 𝜌𝑖,𝑗 𝑡+1       =𝜌𝑖,𝑗𝑡−Δ𝑡3𝜌𝑢𝑖,𝑗 𝑡−4𝜌𝑢𝑖+1,𝑗 𝑡+𝜌𝑢𝑖+2,𝑗 𝑡 2Δ𝑥 ( 3-92) 𝑇𝑎𝑝𝑎 𝑑𝑒𝑟𝑒𝑐ℎ𝑎: 𝜌𝑖,𝑗 𝑡+1       =𝜌𝑖,𝑗𝑡+Δ𝑡3𝜌𝑢𝑖,𝑗−4𝜌𝑢𝑖−1,𝑗+𝜌𝑢𝑖−2,𝑗 2Δ𝑥 ( 3-93) 𝑇𝑎𝑝𝑎 𝑖𝑛𝑓𝑒𝑟𝑖𝑜𝑟: 𝜌𝑖,𝑗 𝑡+1       =𝜌𝑖,𝑗𝑡−Δ𝑡3𝜌𝑣𝑖,𝑗−4𝜌𝑣𝑖,𝑗+1+𝜌𝑣𝑖,𝑗+2 2Δ𝑥 ( 3-94) 𝑇𝑎𝑝𝑎 𝑠𝑢𝑝𝑒𝑟𝑖𝑜𝑟: 𝜌𝑖,𝑗 𝑡+1       =𝜌𝑖,𝑗𝑡+Δ𝑡3𝜌𝑣𝑖,𝑗−4𝜌𝑣𝑖,𝑗−1+𝜌𝑣𝑖,𝑗−2 2Δ𝑥 − 𝜌𝑖+1,𝑗−𝜌𝑖−1,𝑗 2Δ𝑥 ( 3-95) Quedando las densidades en el instante t+1(corrector) como: 𝑇𝑎𝑝𝑎 𝑖𝑧𝑞𝑢𝑖𝑒𝑟𝑑𝑎: 𝜌𝑖,𝑗 𝑡+1=0.5(𝜌𝑖,𝑗𝑡+𝜌𝑖,𝑗𝑡+1       −Δ𝑡3𝜌𝑢𝑖,𝑗 𝑡+1       −4𝜌𝑢𝑖+1,𝑗 𝑡+1       +𝜌𝑢𝑖+2,𝑗 𝑡+1       2Δ𝑥 ) ( 3-96) 𝑇𝑎𝑝𝑎 𝑑𝑒𝑟𝑒𝑐ℎ𝑎: 𝜌𝑖,𝑗 𝑡+1=0.5(𝜌𝑖,𝑗𝑡+𝜌𝑖,𝑗𝑡+1       +Δ𝑡3𝜌𝑢𝑖,𝑗 𝑡+1       −4𝜌𝑢𝑖+1,𝑗 𝑡+1       +𝜌𝑢𝑖+2,𝑗 𝑡+1       2Δ𝑥 ) ( 3-97) 𝑇𝑎𝑝𝑎 𝑖𝑛𝑓𝑒𝑟𝑖𝑜𝑟: 𝜌𝑖,𝑗 𝑡+1=0.5(𝜌𝑖,𝑗𝑡+𝜌𝑖,𝑗𝑡+1       −Δ𝑡3𝜌𝑣𝑖,𝑗 𝑡+1       −4𝜌𝑣𝑖+1,𝑗 𝑡+1       +𝜌𝑣𝑖+2,𝑗 𝑡+1       2Δ𝑥 ) ( 3-98) 38 𝑇𝑎𝑝𝑎 𝑠𝑢𝑝𝑒𝑟𝑖𝑜𝑟: 𝜌𝑖,𝑗 𝑡+1=0.5(𝜌𝑖,𝑗𝑡+𝜌𝑖,𝑗𝑡+1       +Δ𝑡3𝜌𝑣𝑖,𝑗 𝑡+1       −4𝜌𝑣𝑖+1,𝑗 𝑡+1       +𝜌𝑣𝑖+2,𝑗 𝑡+1       2Δ𝑥 − 𝜌𝑖+1,𝑗 𝑡+1       −𝜌𝑖−1,𝑗 𝑡+1       2Δ𝑥 ) ( 3-99) Aplicamos el esquema de MacCormack: 1. Etapa de predicción (Diferencia Progresiva) En esta fase calculamos las variables 𝜌𝑡+1       ,𝜌𝑢𝑡+1       𝑦 𝜌𝑣𝑡+1       mediante las ecuaciones discretizadas de continuidad y de cantidad de movimiento. 2. Decodificación de variables Se calcula u y v a partir de las variables que obtenemos de las ecuaciones. 3. Condiciones de contorno. Ahora aplicaremos las condiciones de contorno en velocidad y en densidad, como se acaba de explicar. 4. Actualización de variables Es importante actualizar nuestras variables auxiliares 𝜌𝑢𝑡+1       ,𝜌𝑣𝑡+1       ,𝜌𝑢2𝑡+1       ,𝜌𝑣2𝑡+1       𝑦 𝜌𝑢𝑣𝑡+1       despues de aplicar las condiciones de contorno. 5. Etapa de corrección (Diferencia regresiva). En esta fase calculamos las variables 𝜌𝑡+1,𝜌𝑢𝑡+1 𝑦 𝜌𝑣𝑡+1 mediante las ecuaciones discretizadas de continuidad y de cantidad de movimiento. 6. Decodificación de variables Se calcula u y v a partir de las variables que obtenemos de las ecuaciones. 7. Condiciones de contorno. Las condiciones de contorno de esta etapa se realizan siguiendo el esquema anterior 8. Actualización de variables Es importante actualizar nuestras variables auxiliares 𝜌𝑢𝑡+1 ,𝜌𝑣𝑡+1,𝜌𝑢2𝑡+1 ,𝜌𝑣2𝑡+1 𝑦 𝜌𝑢𝑣𝑡+1 despues de aplicar las condiciones de contorno. 9. Repetimos el proceso hasta alcanzar el tiempo deseado. Un aspecto que aun no hemos analizado es el incremento de tiempo que vamos a seleccionar. Hemos comentado en el primer apartado que, a pesar de no existir una demostración de una condición de estabilidad para problemas multidimensionales, la condición de Courant-Friedrich-Levy funciona considerablemente bien en la práctica, la recordamos: 𝛥𝑡≤1 |𝑢| 𝛥𝑥+|𝑣| 𝛥𝑦+𝑐√1 𝛥𝑥2+1 𝛥𝑦2 ( 3-100) 39 Incluso mejor, para problemas que conciernen a la resolución de las ecuaciones de Navier-Stokes, Tannehill et al., nos propone la siguiente fórmula empírica. 𝛥𝑡∗≤𝜎Δ𝑡 1+2/𝑅𝑒Δ ( 3-101) Donde 𝑅𝑒Δ es el número de Reynolds basado en el mínimo de Δ𝑥 𝑦 Δ𝑦 y 𝜎=0.9, es un factor de seguridad. Observamos que para valores elevados del Número de Reynolds basado en Δ las condiciones son similares, siendo la segunda algo más restrictiva. Volviendo a la condición CFL observamos que la podemos simplificar si tenemos en cuenta que Δ𝑥 ha sido elegido igual a Δ𝑦 y que el valor máximo de u y v será del mismo orden de magnitud. 𝛥𝑡≤Δ𝑥 2|𝑢|+𝑐√2 ( 3-102) Que observamos que es más exigente que la que obteníamos en una variable. 𝛥𝑡≤Δ𝑥 |𝑢|+𝑐 ( 3-103) A modo de anotación y por si se quisiese optimizar el código R.W. MacCormack hizo una revisión de su método en el que divide su esquema en etapas unidimensionales, lo que como acabamos de ver reduce las restricciones de tiempo y mejora aún más si comtemplamos mallas en las que Δ𝑥 y Δ𝑦 sean muy diferentes ya que se podría avanzar en cada dirección con el incremento de tiempo máximo permitido. No obstante, vamos a aplicar el esquema original por su simplicidad y buenos resultados. Vamos a analizar los casos de 𝑅𝑒=100 𝑦 𝑅𝑒=400 para 𝑀=0.05 para mantenernos siempre en el límite incompresible. 40 Figura( 3-8 ) Figura( 3-9 ) Podemos observar como la tapa ‘arrastra’ al fluido generando la rotación y formando el remolino tan característico. Las líneas de corriente más oscuras corresponden a los resultados de Hou y el campo vectorial de velocidades gris claro al obtenido por Kundu, Cohen & Dowling. 41 Figura( 3-10 ) En el detalle de la imagen se aprecia como el centro de la celda convectiva se encuentra en [0,62 0,74] que son muy próximos a los valores obtenidos por Hou et al.(1995) ,realizados mediante el método latticeBoltzman, que acordaban que el centro se encontraba en [0,6196 0,7373]. Figura( 3-11 ) 48 𝑒𝑐𝑢𝑎𝑐𝑖ó𝑛 𝑑𝑒 𝑒𝑠𝑡𝑎𝑑𝑜: 𝑝= 𝜌 𝛾𝑀2 ( 4-4) En las próximas aplicaciones por simplicidad se escribirán las ecuaciones de Navier-Stokes en notación matricial. El problema está sujeto a las condiciones de contorno: 𝑢(0,𝑦)=1 𝑣(0,𝑦)=0 𝑢(1,𝑦)=0 𝑣(1,𝑦)=0 ( 4-5) 𝑢(𝑥,0)=1 𝑣(𝑥,0)=0 𝑢(𝑥,1)=1 𝑣(𝑥,1)=0 ( 4-6) 𝑢(𝑥𝑐𝑖,𝑦𝑐)=0 𝑣(𝑥𝑐𝑖,𝑦𝑐)=0 𝑢(𝑥𝑐𝑓,𝑦𝑐)=0 𝑣(𝑥𝑐𝑓,𝑦𝑐)=0 ( 4-7) 𝑢(𝑥𝑐,𝑦𝑐𝑖)=0 𝑣(𝑥𝑐,𝑦𝑐𝑖)=0 𝑢(𝑥𝑐,𝑦𝑐𝑓)=0 𝑣(𝑥𝑐,𝑦𝑐𝑓)=0 ( 4-8) Siendo 𝑥𝑐𝑖 la coordenada x del inicio del cuadrado, 𝑥𝑐𝑓 la coordenada x del final del cuadrado y para la coordenada y la nomenclatura es análoga. Es necesario imponer condiciones de contorno en densidad (o presión). Vamos a realizar esto utilizando la ecuación de continuidad: 𝑇𝑎𝑝𝑎 𝑖𝑛𝑓𝑒𝑟𝑖𝑜𝑟(𝑑(𝜌𝑢) 𝑑𝑥 =𝑑𝜌 𝑑𝑥)𝑑𝜌 𝑑𝑡=−𝑑𝜌 𝑑𝑥−𝑑(𝜌𝑣) 𝑑𝑦 ( 4-9) 𝑇𝑎𝑝𝑎 𝑠𝑢𝑝𝑒𝑟𝑖𝑜𝑟(𝑑(𝜌𝑢) 𝑑𝑥 =𝑑𝜌 𝑑𝑥)𝑑𝜌 𝑑𝑡=−𝑑𝜌 𝑑𝑥−𝑑(𝜌𝑣) 𝑑𝑦 ( 4-10) 𝑇𝑎𝑝𝑎 𝑖𝑧𝑞𝑢𝑖𝑒𝑟𝑑𝑎(𝑑(𝜌𝑣) 𝑑𝑦 =0)𝑑𝜌 𝑑𝑡=−𝑑(𝜌𝑢) 𝑑𝑥 ( 4-11) 𝑇𝑎𝑝𝑎 𝑑𝑒𝑟𝑒𝑐ℎ𝑎(𝑑(𝜌𝑣) 𝑑𝑦 =0)𝑑𝜌 𝑑𝑡=−𝑑(𝜌𝑢) 𝑑𝑥 ( 4-12) Que observamos que si aplicamos las diferencias progresivas, regresivas y centrales obtenidas en el apartado 3.2 podremos calcular las densidades en la frontera de nuestro canal. Para las condiciones de densidad en el cuadrado, se puede actualizar o bien con la ecuación de continuidad: 𝐿𝑎𝑑𝑜 𝑖𝑛𝑓𝑒𝑟𝑖𝑜𝑟(𝑑(𝜌𝑢) 𝑑𝑥 =0)𝑑𝜌 𝑑𝑡=−𝑑(𝜌𝑣) 𝑑𝑦 ( 4-13) 𝐿𝑎𝑑𝑜 𝑠𝑢𝑝𝑒𝑟𝑖𝑜𝑟(𝑑(𝜌𝑢) 𝑑𝑥 =0)𝑑𝜌 𝑑𝑡=−𝑑(𝜌𝑣) 𝑑𝑦 ( 4-14) 49 𝐿𝑎𝑑𝑜 𝑖𝑧𝑞𝑢𝑖𝑒𝑟𝑑𝑜(𝑑(𝜌𝑣) 𝑑𝑦 =0)𝑑𝜌 𝑑𝑡=−𝑑(𝜌𝑢) 𝑑𝑥 ( 4-15) 𝑇𝑎𝑝𝑎 𝑑𝑒𝑟𝑒𝑐ℎ𝑜(𝑑(𝜌𝑣) 𝑑𝑦 =0)𝑑𝜌 𝑑𝑡=−𝑑(𝜌𝑢) 𝑑𝑥 ( 4-16) O bien desde la ecuación de cantidad de movimiento, que es más compleja, pero en este caso parece dar mejores resultados: 𝐿𝑎𝑑𝑜 𝑖𝑛𝑓𝑒𝑟𝑖𝑜𝑟 𝑑𝜌 𝑑𝑦=𝑀2(1 𝑅𝑒(43𝑑2𝑣 𝑑𝑦2+13𝑑2𝑢 𝑑𝑥𝑑𝑦+𝑑2𝑣 𝑑𝑥2)−𝑑(𝜌𝑢𝑣) 𝑑𝑥 −𝑑(𝜌𝑣2) 𝑑𝑦 −𝑑(𝜌𝑣) 𝑑𝑡 ) =𝑀2 𝑅𝑒(43𝑑2𝑣 𝑑𝑦2+13𝑑2𝑢 𝑑𝑥𝑑𝑦) ( 4-17) 𝐿𝑎𝑑𝑜 𝑠𝑢𝑝𝑒𝑟𝑖𝑜𝑟 𝑑𝜌 𝑑𝑦=𝑀2(1 𝑅𝑒(43𝑑2𝑣 𝑑𝑦2+13𝑑2𝑢 𝑑𝑥𝑑𝑦+𝑑2𝑣 𝑑𝑥2)−𝑑(𝜌𝑢𝑣) 𝑑𝑥 −𝑑(𝜌𝑣2) 𝑑𝑦 −𝑑(𝜌𝑣) 𝑑𝑡 ) =𝑀2 𝑅𝑒(43𝑑2𝑣 𝑑𝑦2+13𝑑2𝑢 𝑑𝑥𝑑𝑦) ( 4-18) 𝐿𝑎𝑑𝑜 𝑖𝑧𝑞𝑢𝑖𝑒𝑟𝑑𝑜 𝑑𝜌 𝑑𝑥=𝑀2(1 𝑅𝑒(43𝑑2𝑢 𝑑𝑥2+13𝑑2𝑣 𝑑𝑥𝑑𝑦+𝑑2𝑢 𝑑𝑦2)−𝑑(𝜌𝑢2) 𝑑𝑥 −𝑑(𝜌𝑣𝑢) 𝑑𝑦 −𝑑(𝜌𝑢) 𝑑𝑡 ) =𝑀2 𝑅𝑒(43𝑑2𝑢 𝑑𝑥2+13𝑑2𝑣 𝑑𝑥𝑑𝑦) ( 4-19) 𝐿𝑎𝑑𝑜 𝑑𝑒𝑟𝑒𝑐ℎ𝑜 𝑑𝜌 𝑑𝑥=𝑀2(1 𝑅𝑒(43𝑑2𝑢 𝑑𝑥2+13𝑑2𝑣 𝑑𝑥𝑑𝑦+𝑑2𝑢 𝑑𝑦2)−𝑑(𝜌𝑢2) 𝑑𝑥 −𝑑(𝜌𝑣𝑢) 𝑑𝑦 −𝑑(𝜌𝑢) 𝑑𝑡 ) =𝑀2 𝑅𝑒(43𝑑2𝑢 𝑑𝑥2+13𝑑2𝑣 𝑑𝑥𝑑𝑦) ( 4-20) Aplicando las aproximaciones establecidas en el capítulo 3 se pueden despejar las variables en el contorno del cuadrado. Para calcular los coeficientes de resistencia y sustentación vamos a utilizar las siguientes expresiones: 𝐶𝑙=𝑙 12𝜌𝑈2𝐷=∫(𝑃𝑖𝑛𝑓−𝑃𝑠𝑢𝑝)𝑑𝑥+∫(𝑓𝑦𝑖𝑧𝑞+𝑓𝑦𝑑𝑒𝑟)𝑑𝑦+∫(𝑓𝑦𝑠𝑢𝑝+𝑓𝑦𝑖𝑛𝑓)𝑑𝑥 𝑖2 𝑖1 𝑗2 𝑗1 𝑖2 𝑖112𝜌𝑈2𝐷 ( 4-21) 50 Siendo 𝑓𝑦 los esfuerzos tangenciales que obtenemos en las paredes: 𝑓𝑦𝑖𝑧𝑞=−𝜇𝑑𝑣 𝑑𝑥𝑥=𝑖1 ( 4-22) 𝑓𝑦𝑑𝑒𝑟=𝜇𝑑𝑣 𝑑𝑥𝑥=𝑖2 ( 4-23) 𝑓𝑦𝑑𝑒𝑟=𝜇𝑑𝑣 𝑑𝑥𝑥=𝑖2 ( 4-24) 𝑓𝑦𝑠𝑢𝑝=2𝜇𝑑𝑣 𝑑𝑦𝑦=𝑗2 ( 4-25) 𝑓𝑦𝑖𝑛𝑓=−2𝜇𝑑𝑣 𝑑𝑦𝑦=𝑗1 ( 4-26) Convirtiendo a nuestras variables adimensionales: 𝐶𝑙 =𝜌𝑈2𝐷∫(𝑃𝑖𝑛𝑓 ∗−𝑃𝑠𝑢𝑝 ∗)𝑑𝑥∗+𝜇𝑈∫(𝑓𝑦∗𝑖𝑧𝑞+𝑓𝑦∗𝑑𝑒𝑟)𝑑𝑦∗ 𝑗2 𝐷 𝑗1 𝐷+𝜇𝑈∫(𝑓𝑦∗𝑠𝑢𝑝+𝑓𝑦∗𝑖𝑛𝑓)𝑑𝑥∗ 𝑗2 𝐷 𝑗1 𝐷 𝑖2 𝐷 𝑖1 𝐷12𝜌𝑈2𝐷 =2∫(𝑃𝑖𝑛𝑓 ∗−𝑃𝑠𝑢𝑝 ∗)𝑑𝑥∗+2 𝑅𝑒∫ (𝑓𝑦∗𝑖𝑧𝑞+𝑓𝑦∗𝑑𝑒𝑟)𝑑𝑦∗+2 𝑅𝑒∫ (𝑓𝑦∗𝑠𝑢𝑝+𝑓𝑦∗𝑖𝑛𝑓)𝑑𝑥∗ 𝑖2 𝐷 𝑖1 𝐷 𝑗2 𝐷 𝑗1 𝐷 𝑖2 𝐷 𝑖1 𝐷 ( 4-27) Lo cual tiene sentido físico ya que un aumento de la diferencia de presiones entre las caras superior e inferior contribuye positivamente a la sustentación. También es coherente que aumentar el número de Reynolds produzca una perdida de importancia de las fuerzas de viscosidad. Si tenemos en cuenta que 𝑃∗=𝜌∗ 𝑀2: 𝐶𝑙=2 𝑀2∫(𝜌𝑖𝑛𝑓 ∗−𝜌𝑠𝑢𝑝 ∗)𝑑𝑥∗+ 𝑖2 𝐷 𝑖1 𝐷2 𝑅𝑒∫ (𝑓𝑦∗𝑖𝑧𝑞+𝑓𝑦∗𝑑𝑒𝑟)𝑑𝑦∗+2 𝑅𝑒∫ (𝑓𝑦∗𝑠𝑢𝑝+𝑓𝑦∗𝑖𝑛𝑓)𝑑𝑥∗ 𝑖2 𝐷 𝑖1 𝐷 𝑗2 𝐷 𝑗1 𝐷 ( 4-28) Análogamente para el coeficiente de resistencia: 51 𝐶𝑑=𝑑 12𝜌𝑈2𝐷=∫(𝑃𝑖𝑧𝑞−𝑃𝑑𝑒𝑟)𝑑𝑦+∫(𝑓𝑥𝑖𝑧𝑞+𝑓𝑥𝑑𝑒𝑟)𝑑𝑦+∫(𝑓𝑥𝑠𝑢𝑝+𝑓𝑥𝑖𝑛𝑓)𝑑𝑥 𝑖2 𝑖1 𝑗2 𝑗1 𝑗2 𝑗112𝜌𝑈2𝐷 ( 4-29) Siendo 𝑓𝑥 los esfuerzos tangenciales que obtenemos en las paredes: 𝑓𝑥𝑖𝑧𝑞=−2𝜇𝑑𝑢 𝑑𝑥𝑥=𝑖1 ( 4-30) 𝑓𝑥𝑑𝑒𝑟=2𝜇𝑑𝑢 𝑑𝑥𝑥=𝑖2 ( 4-31) 𝑓𝑥𝑠𝑢𝑝=𝜇𝑑𝑢 𝑑𝑦𝑦=𝑗2 ( 4-32) 𝑓𝑥𝑖𝑛𝑓=−𝜇𝑑𝑢 𝑑𝑦𝑦=𝑗1 ( 4-33) Convirtiendo a nuestras variables adimensionales: 𝐶𝑑 =𝜌𝑈2𝐷∫(𝑃𝑖𝑧𝑞 ∗−𝑃𝑑𝑒𝑟 ∗)𝑑𝑦∗+𝜇𝑈∫(𝑓𝑥∗𝑖𝑧𝑞+𝑓𝑥∗𝑑𝑒𝑟)𝑑𝑦∗ 𝑗2 𝐷 𝑗1 𝐷+𝜇𝑈∫(𝑓𝑥∗𝑠𝑢𝑝+𝑓𝑥∗𝑖𝑛𝑓)𝑑𝑥∗ 𝑗2 𝐷 𝑗1 𝐷 𝑗2 𝐷 𝑗1 𝐷12𝜌𝑈2𝐷 =2∫ (𝑃𝑖𝑧𝑞 ∗−𝑃𝑑𝑒𝑟 ∗)𝑑𝑦∗+2 𝑅𝑒∫ (𝑓𝑥∗𝑖𝑧𝑞+𝑓𝑥∗𝑑𝑒𝑟)𝑑𝑦∗+2 𝑅𝑒∫ (𝑓𝑥∗𝑠𝑢𝑝+𝑓𝑥∗𝑖𝑛𝑓)𝑑𝑥∗ 𝑖2 𝐷 𝑖1 𝐷 𝑗2 𝐷 𝑗1 𝐷 𝑗2 𝐷 𝑗1 𝐷 ( 4-34) Lo cual tiene sentido físico ya que un aumento de la diferencia de presiones entre las caras izquierda y derecha contribuye positivamente a la resistencia. También es coherente que aumentar el número de Reynolds produzca una perdida de importancia de las fuerzas de viscosidad. Si tenemos en cuenta que 𝑃∗=𝜌∗ 𝑀2: 𝐶𝑑=2 𝑀2∫(𝜌𝑖𝑧𝑞 ∗−𝜌𝑑𝑒𝑟 ∗)𝑑𝑦∗+ 𝑗2 𝐷 𝑗1 𝐷2 𝑅𝑒∫ (𝑓𝑥∗𝑖𝑧𝑞+𝑓𝑥∗𝑑𝑒𝑟)𝑑𝑦∗+2 𝑅𝑒∫ (𝑓𝑥∗𝑠𝑢𝑝+𝑓𝑥∗𝑖𝑛𝑓)𝑑𝑥∗ 𝑖2 𝐷 𝑖1 𝐷 𝑗2 𝐷 𝑗1 𝐷 ( 4-35) Comenzaremos analizando el flujo alrededor del cuadrado para 𝑅𝑒=20 𝑦 𝑀=0.05 Antes de analizar los coeficientes de resistencia y sustentación vamos a analizar el patrón de flujo que se obtiene entorno a la caja y que le ocurre a la corriente al rebordear a esta. 52 Figura( 4-2 ) Figura( 4-3 ) Observamos como inmediatamente después de bordear la caja se forman dos torbellinos que permanecen adheridos, y observamos como en ningún caso se produce oscilación en la estela, obteniendo un régimen estacionario. Vamos a representar ahora los coeficientes de resistencia y sustentación. 53 Figura( 4-4 ) Observamos como el valor no es idénticamente nulo, esto se debe a asímetrias en la malla, no obstante, se puede observar que si refinamos la malla este valor va disminuyendo hasta aproximarse mucho a 0. Figura( 4-5 ) Percibimos con esta imagen como estas asimetrías se han reducido considerablemente hasta un valor de 0.03, ya muy próximo a sustentación nula que era el comportamiento esperado. 54 Figura( 4-6 ) En la figura 4-6 observamos como el valor del Cd no es residual como ocurria con el Cl en este caso como era de esperar existe una resistencia importante al movimiento ya que la corriente se desprende, provocando las burbujas de recirculación en sotavento produciendo una importante resistencia. Otro aspecto significativo es el hecho de que para este número de Reynolds aparezca un régimen estacionario, algo que no será verdad para cualquier número de Reynolds como veremos a continuación. Re=100. Figura( 4-7 ) 55 Figura( 4-8 ) Figura( 4-9 ) En la secuencia de imágenes anterior,ordenadas temporalmente, se muestra como el fluido comienza a rotar al rebordear la caja formándose un torbellino que se va convectando aguas abajo haciendo que la estela oscile lo que como veremos a continuación provocará variaciones temporales en los coeficientes de resistencia y de sustentación. 56 Figura( 4-10 ) Como comentábamos antes, obtenemos una oscilación en el coeficiente de sustentación, esto se debe a la oscilación de la estela como se ha mostrado anteriormente, aunque si que observamos que debido a la simetría del problema esta oscilación es simétrica respecto a la horizontal de una amplitud aproximadamente 0.8, muy similar a la obtenida por Kundu, Cohen &Dowling, que obtenían un valor de 0.77. Figura( 4-11 ) 57 Observamos que como era de esperar la oscilación de la estela produce una variación de la resistencia, no obstante, esta es mucho más reducida en amplitud y de frecuencia mucho más alta, observamos como su valor medio este entorno a los 3.4, en la siguiente imagen se aprecia con más claridad. De nuevo un valor muy similar al predicho por Kundu, Cohen &Dowling que obtenían un valor medio de 3.35. Figura( 4-12 ) Es interesante comparar este valor medio obtenido con el del apartado anterior y observar que se ha reducido prácticamente a la mitad, lo cual tiene sentido físico ya que al aumentar el número de Reynolds disminuye la importancia relativa de las fuerzas de viscosidad y por tanto de la resistencia. Vamos a realizar un análisis de sensibilidad de la relación de aspecto de la cavidad respecto del cuadrado para un valor del número de Reynolds de 20 y después de un tiempo adimensional de 15, cuando ya se ha alcanzado el régimen estacionario, en relación al coeficiente de resistencia. La posición relativa del cuadrado se ha mantenido constante. . Tabla ( 4-1 ) 𝑳/𝑫 𝑯/𝑫 𝑪𝒅 35 3 5.5703 100 3 6.0650 La condición de contorno que se ha puesto en la salida del conducto es de velocidad uniforme igual a velocidad de entrada, por tanto, al aumentar la longitud de la cavidad se tiende al caso real ya que los efectos perturbadores producto de rebordear el objeto ya son despreciables. En cambio si reducimos la longitud de la cavidad obligamos a los efectos perturbadores de la estela a amortiguarse antes de lo que lo harían en un caso real. 64 Figura( 4-22 ) Observamos como en la figura 4-15 correspondiente a 𝑅𝑒=2 la corriente permanece adherida, mientras que en la figura 4-16 correspondiente a 𝑅𝑒=5 esto ya no ocurre y se comienzan a formar torbellinos que se mantienen adheridos a sotavento, sin provocar la oscilación de la estela. En la figura 4-17 ya se observan los torbellinos de mayor intensidad pero se mantienen adheridos. En la figura 4-18 ya observamos como se comienza a romper la simetría y se comienza a emitir vorticidad aguas abajo, el valor del número de Reynolds crítico para el cual comienza la emisión de torbellinos será muy próximo pero ligeramente inferiror a 𝑅𝑒=40 . Finalmente en la figura 4-19 se aprecia muy claramente como se están emitiendo torbellinos aguas abajo provocando la oscilación de la estela que antes comentábamos. Observamos como los resultados son muy parejos a los del cilindro circular infinito, probablemente en nuestro caso los umbrales sean ligeramente inferiores, debido a que los angulos rectos en los que se unen las caras de nuestro cilindro facilitan el desprendimiento de la corriente. Llegados a este punto nos surge la siguiente preocupación, ¿Con que frecuencia se desprenden torbellinos para distintos números de Reynolds? Como hemos visto anteriormente esto se puede estudiar mediante el número de Strouhal que es equivalente a una frecuencia adimensional. La adimensionalización que hemos realizado no puede ser más adecuada para el estudio de este fenómeno ya que para calcular la frecuencia adimensional y por ende el número de Strouhal únicamente debemos calcular el periodo adimensional he invertirlo, esto se puede hacer por ejemplo observando las oscilaciones del coeficiente de sustentación y obteniendo su periodo. En la figura 4-22 observamos los resultados obtenidos mediante el esquema de MacCormack para una cavidad de relaciones de aspecto 𝐿𝐷=35 y 𝐻 𝐷=6 los cuales deberían aproximar considerablemente bien a un cilindro inmerso en una corriente uniforme junto con los resultados de la bibliografía. 65 Figura( 4-23 ) Donde se observa como el acuerdo con los datos obtenidos de diverosos autores es muy elevado. Para obtener una precisión incluso mayor podíamos ampliar el dominio computacional e imponer las condiciones de contorno más lejos del cilindro. Hemos optado por esta solución de compromiso porque proporcionaba resultados que se ajustaban a la bibliografía sin requerir un tiempo de computación excesivamente elevado. Como hemos comentado en la introducción en la ingeniería la variación de sustentación puede producir problemas estructurales en antenas de vehículos, torres de refrigeración, rascacielos, en definitiva estructuras que vean una corriente aproximadamente uniforme. Para evitar la rotura por fatiga de estos elementos suelen incluir aletas para tratar de reducir estas oscilaciones ya que si se evite que los torbellinos estén en contacto directo la inestabilidad se reduce considerablemente. A continuación mostramos dos imágenes una correspondiente al flujo alrededor de un cilindro circular, y en la siguiente el cilindro incorpora una aleta para evitar la formación de vórtices. Figura( 4-24 ) Se aprecia claramente la reducción en la oscilación de la estela. Vamos a mostrar el efecto para nuestro cuadrado. 66 Figura( 4-25 ) Figura( 4-26 ) 67 Figura( 4-27 ) Se aprecia claramente como la inclusión de la aleta proporciona una estela mucho menos oscilante. En la aleta mas corta observamos como ambos vórtices se encuentran generando la oscilación de la estela pero con menor intensidad que en el caso sin aleta. En el caso de la aleta mas larga, figura 4-26 observamos como los torbellinos apenas están en contacto generando una oscilación mínima de la estela. En la siguiente página mostramos el código MATLAB utilizado para la resolución de estos problemas. 68 %=====================FLUJO ALREDEDOR DE UN PRISMA DE SECCIÓN CUADRADA=================% %================================PARÁMETROS================================= ===========% close all;clear all; nx=349;%intervalos de partición en x ny=59;%intervalos de partición en y nt=20000;%intervalos de partición en t D=1; %ancho OBSTACULO Dt=0.002; %incremento temporal Dx=35*D/(nx); %incremento eje x Dy=6*D/(ny); %incremento eje y U=1; %velocidad cavidad %===================================================================% %===============NODOS DEL OBSTACULO=================================% i1=round(15/35*(nx+1)); i2=round(16/35*(nx+1)); j1=round(2.5/6*(ny+1)); j2=round(3.5/6*(ny+1)); %===================================================================% %===============NUMEROS ADIMENSIONALES==============================% Re=40; %numero de Reynolds M=0.05; %numero de Mach %===================COEFICIENTES====================================% a1=Dt/Dx; a2=Dt/Dy; a3=Dt/(Dx*M^2); a4=Dt/(Dy*M^2); a5=4*Dt/(3*Re*Dx^2); a6=Dt/(Re*Dy^2); a7=Dt/(Re*Dx^2); a8=4*Dt/(3*Re*Dy^2); a9=Dt/(12*Re*Dx*Dy); a10=2*(a5+a6); a11=2*(a7+a8); a12=8/9*M^2/Dx/Re; a13=M^2/18/Dy/Re; %=================CONDICIONES INICIALES=============================% P_u=zeros(nx+1,ny+1);%velocidad x en los puntos del mallado P_v=zeros(nx+1,ny+1);%velocidad y en los puntos del mallado P_rho=ones(nx+1,ny+1);% densidad en los puntos del mallado P_rhou=zeros(nx+1,ny+1);% densidad*u en los puntos del mallado P_rhov=zeros(nx+1,ny+1);%densidad*v en los puntos del mallado P_rhou2=zeros(nx+1,ny+1);% densidad*u^2 en los puntos del mallado P_rhov2=zeros(nx+1,ny+1);%densidad*v^2 en los puntos del mallado P_rhouv=zeros(nx+1,ny+1);%densidad*uv en los puntos del mallado u=zeros(nx+1,ny+1);%velocidad x en los puntos del mallado v=zeros(nx+1,ny+1);%velocidad y en los puntos del mallado rho=ones(nx+1,ny+1);% densidad en los puntos del mallado rhou=zeros(nx+1,ny+1);% densidad*u en los puntos del mallado rhov=zeros(nx+1,ny+1);%densidad*v en los puntos del mallado rhou2=rhou.*u(:,:);% densidad*u^2 en los puntos del mallado rhov2=zeros(nx+1,ny+1);%densidad*v^2 en los puntos del mallado rhouv=zeros(nx+1,ny+1);%densidad*uv en los puntos del mallado %========================BUCLE TEMPORAL===============================% cont=1; k=0; for q=1:nt 69 k=k+1; q %=========================PREDICTOR==================================% P_rho(2:nx,2:ny)=rho(2:nx,2:ny)-a1*(rhou(3:nx+1,2:ny)-rhou(2:nx,2:ny))- a2*(rhov(2:nx,3:ny+1)-rhov(2:nx,2:ny)); P_rhou(2:nx,2:ny)=rhou(2:nx,2:ny)-a3*(rho(3:nx+1,2:ny)-rho(2:nx,2:ny))- a1*(rhou2(3:nx+1,2:ny)-rhou2(2:nx,2:ny))-a2*(rhouv(2:nx,3:ny+1)- rhouv(2:nx,2:ny))-a10*u(2:nx,2:ny)+a5*(u(3:nx+1,2:ny)+u(1:nx1,2:ny))+a6*(u(2:nx,3:ny+1)+u(2:nx,1:ny-1))+a9*(v(3:nx+1,3:ny+1)+v(1:nx1,1:ny-1)-v(1:nx-1,3:ny+1)-v(3:nx+1,1:ny-1)); P_rhov(2:nx,2:ny)=rhov(2:nx,2:ny)-a4*(rho(2:nx,3:ny+1)-rho(2:nx,2:ny))- a1*((rhouv(3:nx+1,2:ny))-(rhouv(2:nx,2:ny)))-a2*(rhov2(2:nx,3:ny+1)- rhov2(2:nx,2:ny))-a11*v(2:nx,2:ny)+a7*(v(3:nx+1,2:ny)+v(1:nx1,2:ny))+a8*(v(2:nx,3:ny+1)+v(2:nx,1:ny-1))+a9*(u(3:nx+1,3:ny+1)+u(1:nx1,1:ny-1)-u(1:nx-1,3:ny+1)-u(3:nx+1,1:ny-1)); %=========================DECODIFICACIÓN=============================% P_u(2:nx,2:ny)=P_rhou(2:nx,2:ny)./P_rho(2:nx,2:ny); P_v(2:nx,2:ny)=P_rhov(2:nx,2:ny)./P_rho(2:nx,2:ny); %=============================C.C====================================% P_rho(1,1:ny+1)=rho(1,1:ny+1)-(a1/2)*(-rhou(3,1:ny+1)+4*rhou(2,1:ny+1)- 3*rhou(1,1:ny+1));%%%%x=0 densidad P_rho(nx+1,1:ny+1)=rho(nx+1,1:ny+1)+(a1/2)*(-rhou(nx1,1:ny+1)+4*rhou(nx,1:ny+1)-3*rhou(nx+1,1:ny+1));%%%x=d densidad P_rho(2:nx,1)=rho(2:nx,1)-(a2/2)*(-rhov(2:nx,3)+4*rhov(2:nx,2)- 3*rhov(2:nx,1))-(a1*U/2)*(rho(3:nx+1,1)-rho(1:nx-1,1));%%%%y=0 densidad P_rho(2:nx,ny+1)=rho(2:nx,ny+1)+a2/2*(-rhov(2:nx,ny-1)+4*rhov(2:nx,ny)- 3*rhov(2:nx,ny+1))-(a1*U/2)*(rho(3:nx+1,ny+1)-rho(1:nx-1,ny+1));%%%y=d densidad P_rho(i1,j1:j2)=1/3*(4*rho(i1-1,j1:j2)-rho(i1-2,j1:j2))+a12*(-5*u(i11,j1:j2)+4*u(i1-2,j1:j2)-u(i1-3,j1:j2))-a13*(-(v(i1-2,j1+1:j2+1)-v(i1,j11:j2-1))+4*(v(i1-1,j1+1:j2+1)-v(i1-1,j1-1:j2-1))-3*(v(i1,j1+1:j2+1)-v(i1,j11:j2-1)));%%%%x=0 densidad cubo P_rho(i2,j1:j2)=1/3*(4*rho(i2+1,j1:j2)-rho(i2+2,j1:j2))-a12*(- 5*u(i2+1,j1:j2)+4*u(i2+2,j1:j2)-u(i2+3,j1:j2))-a13*(-(v(i2+2,j1+1:j2+1)- v(i2,j1-1:j2-1))+4*(v(i2+1,j1+1:j2+1)-v(i2+1,j1-1:j2-1))-3*(v(i2,j1+1:j2+1)- v(i2,j1-1:j2-1)));%%%x=d densidad cubo P_rho(i1:i2,j1)=1/3*(4*rho(i1:i2,j1-1)-rho(i1:i2,j1-2))+a12*(- 5*v(i1:i2,j1-1)+4*v(i1:i2,j1-2)-v(i1:i2,j1-3))-a13*(-(u(i1+1:i2+1,j1-2)- u(i1-1:i2-1,j1-2))+4*(u(i1+1:i2+1,j1-1)-u(i1-1:i2-1,j1-1))- 3*(u(i1+1:i2+1,j1)-u(i1-1:i2-1,j1)));%%%%y=0 densidad cubo P_rho(i1:i2,j2)=1/3*(4*rho(i1:i2,j2+1)-rho(i1:i2,j2+2))-a12*(- 5*v(i1:i2,j2+1)+4*v(i1:i2,j2+2)-v(i1:i2,j2+3))-a13*(-(u(i1+1:i2+1,j2+2)- u(i1-1:i2-1,j2+2))+4*(u(i1+1:i2+1,j2+1)-u(i1-1:i2-1,j2+1))- 3*(u(i1+1:i2+1,j2)-u(i1-1:i2-1,j2)));%%%y=d densidad cubo 70 P_u(2:nx,ny+1)=U;%%%%velocidad u y=d P_v(2:nx,ny+1)=0;%%%%velocidad v y=d P_u(2:nx,1)=U;%%%%velocidad u y=0 P_v(2:nx,1)=0;%%%%velocidad v y=0 P_u(nx+1,:)=U;%%%%velocidad u x=d P_v(nx+1,:)=0;%%%%velocidad v x=d P_u(1,:)=U;%%%%velocidad u x=0 P_v(1,:)=0;%%%%velocidad v x=0 P_u(i1:i2,j1:j2)=0; P_v(i1:i2,j1:j2)=0; %===========================ACTUALIZACIÓN========================% P_rhou=P_rho.*P_u; P_rhov=P_rho.*P_v; P_rhou2=P_rhou.*P_u; P_rhouv=P_rhou.*P_v; P_rhov2=P_rhov.*P_v; %==========================CORRECTOR==============================% rho(2:nx,2:ny)=(rho(2:nx,2:ny)+P_rho(2:nx,2:ny) -a1*(P_rhou(2:nx,2:ny)- P_rhou(1:nx-1,2:ny))-a2*(P_rhov(2:nx,2:ny)-P_rhov(2:nx,1:ny-1)))*0.5; rhou(2:nx,2:ny)=0.5*(rhou(2:nx,2:ny)+P_rhou(2:nx,2:ny)- a3*(P_rho(2:nx,2:ny)-P_rho(1:nx-1,2:ny))-a1*(P_rhou2(2:nx,2:ny)- P_rhou2(1:nx-1,2:ny))-a2*(P_rhouv(2:nx,2:ny)-P_rhouv(2:nx,1:ny-1))- a10*P_u(2:nx,2:ny)+a5*(P_u(3:nx+1,2:ny)+P_u(1:nx1,2:ny))+a6*(P_u(2:nx,3:ny+1)+P_u(2:nx,1:ny1))+a9*(P_v(3:nx+1,3:ny+1)+P_v(1:nx-1,1:ny-1)-P_v(1:nx-1,3:ny+1)- P_v(3:nx+1,1:ny-1))); rhov(2:nx,2:ny)= 0.5*(rhov(2:nx,2:ny)+P_rhov(2:nx,2:ny)- a4*(P_rho(2:nx,2:ny)-P_rho(2:nx,1:ny-1))-a1*(P_rhouv(2:nx,2:ny)- P_rhouv(1:nx-1,2:ny))-a2*(P_rhov2(2:nx,2:ny)-P_rhov2(2:nx,1:ny-1))- a11*P_v(2:nx,2:ny)+a7*(P_v(3:nx+1,2:ny)+P_v(1:nx1,2:ny))+a8*(P_v(2:nx,3:ny+1)+P_v(2:nx,1:ny1))+a9*(P_u(3:nx+1,3:ny+1)+P_u(1:nx-1,1:ny-1)-P_u(1:nx-1,3:ny+1)- P_u(3:nx+1,1:ny-1))); %=========================DECODIFICACIÓN=============================% u=rhou./rho; v=rhov./rho; %=============================C.C====================================% rho(1,2:ny)=(P_rho(1,2:ny)+rho(1,2:ny)-a1/2*(- P_rhou(3,2:ny)+4*P_rhou(2,2:ny)-3*P_rhou(1,2:ny)))*0.5;%%%%x=0 rho(nx+1,2:ny)=(P_rho(nx+1,2:ny)+rho(1,2:ny)+a1/2*(-P_rhou(nx1,2:ny)+4*P_rhou(nx,2:ny)-3*P_rhou(nx+1,2:ny)))*0.5;%%%x=d rho(2:nx,1)=(P_rho(2:nx,1)+rho(2:nx,1)-a2/2*(- P_rhov(2:nx,3)+4*P_rhov(2:nx,2)-3*P_rhov(2:nx,1))-a1*U/2*(P_rho(3:nx+1,1)- P_rho(1:nx-1,1)))*0.5;%%%%y=0 rho(2:nx,1+ny)=(P_rho(2:nx,1+ny)+rho(2:nx,1+ny)+a2/2*(-P_rhov(2:nx,ny1)+4*P_rhov(2:nx,ny)-3*P_rhov(2:nx,ny+1))-a1*U/2*(P_rho(3:nx+1,ny+1)- P_rho(1:nx-1,ny+1)))*0.5;%%%y=d 71 rho(i1,j1:j2)=(rho(i1,j1:j2)+1/3*(4*P_rho(i1-1,j1:j2)-P_rho(i12,j1:j2))+a12*(-5*P_u(i1-1,j1:j2)+4*P_u(i1-2,j1:j2)-P_u(i1-3,j1:j2))-a13*(- (P_v(i1-2,j1+1:j2+1)-P_v(i1,j1-1:j2-1))+4*(P_v(i1-1,j1+1:j2+1)-P_v(i1-1,j11:j2-1))-3*(P_v(i1,j1+1:j2+1)-P_v(i1,j1-1:j2-1))))*0.5;%%%%x=0 densidad cubo rho(i2,j1:j2)=(rho(i2,j1:j2)+1/3*(4*P_rho(i2+1,j1:j2)- P_rho(i2+2,j1:j2))-a12*(-5*P_u(i2+1,j1:j2)+4*P_u(i2+2,j1:j2)- P_u(i2+3,j1:j2))-a13*(-(P_v(i2+2,j1+1:j2+1)-P_v(i2,j1-1:j21))+4*(P_v(i2+1,j1+1:j2+1)-P_v(i2+1,j1-1:j2-1))-3*(P_v(i2,j1+1:j2+1)- P_v(i2,j1-1:j2-1))))*0.5;%%%x=d densidad cubo rho(i1:i2,j1)=(rho(i1:i2,j1)+1/3*(4*P_rho(i1:i2,j1-1)-P_rho(i1:i2,j12))+a12*(-5*P_v(i1:i2,j1-1)+4*P_v(i1:i2,j1-2)-P_v(i1:i2,j1-3))-a13*(- (P_u(i1+1:i2+1,j1-2)-P_u(i1-1:i2-1,j1-2))+4*(P_u(i1+1:i2+1,j1-1)-P_u(i11:i2-1,j1-1))-3*(P_u(i1+1:i2+1,j1)-P_u(i1-1:i2-1,j1))))*0.5;%%%%y=0 densidad cubo rho(i1:i2,j2)=(rho(i1:i2,j2)+1/3*(4*P_rho(i1:i2,j2+1)- P_rho(i1:i2,j2+2))-a12*(-5*P_v(i1:i2,j2+1)+4*P_v(i1:i2,j2+2)- P_v(i1:i2,j2+3))-a13*(-(P_u(i1+1:i2+1,j2+2)-P_u(i1-1:i21,j2+2))+4*(P_u(i1+1:i2+1,j2+1)-P_u(i1-1:i2-1,j2+1))-3*(P_u(i1+1:i2+1,j2)- P_u(i1-1:i2-1,j2))))*0.5;%%%y=d densidad cubo u(2:nx,ny+1)=U;%%%%velocidad u y=d v(2:nx,ny+1)=0;%%%%velocidad v y=d u(2:nx,1)=U;%%%%velocidad u y=0 v(2:nx,1)=0;%%%%velocidad v y=0 u(nx+1,:)=U;%%%%velocidad u x=d v(nx+1,:)=0;%%%%velocidad v x=d u(1,:)=U;%%%%velocidad u x=0 v(1,:)=0;%%%%velocidad v x=0 u(i1:i2,j1:j2)=0; v(i1:i2,j1:j2)=0; %===========================ACTUALIZACIÓN========================% rhou=rho.*u; rhov=rho.*v; rhou2=rhou.*u; rhouv=rhou.*v; rhov2=rhov.*v; %======================REPRESENTACIÓN GRÁFICA====================% if k==100 pel(:,:,cont)=[u v]; k=0; cont=cont+1; end %======================CÁLCULO DE COEFICIENTES===================% CL(q)=2*sum(rho(i1:i2,j1)-rho(i1:i2,j2))*Dx/M^2; CD(q)=2*sum(rho(i1,j1:j2)-rho(i2,j1:j2))*Dy/M^2; 72 end figure T=sqrt(v'.^2+u'.^2); quiver(linspace(0,35,nx+1),linspace(0,6,ny+1),u'./T,v'./T) t=Dt*nt 73 5 APLICACIÓN: FLUJO BIDIMENSIONAL ALREDEDOR DE UN CUADRADO CON TEMPERATURA El esquema de esta aplicación es el mismo que la anterior, un cuadrado desplazándose hacia la izquierda en un canal, pero en este caso la temperatura del cuadrado será mayor a la temperatura ambiente permitiéndonos asi estudiar las capas limite Térmicas que pudiesen generarse. Figura( 5-1 ) Las ecuaciones en esta aplicación lógicamente deben de ser complementadas con la ecuación de la energía y condicionesde contorno en temperatura. A continuación se muestran las ecuaciones completas que expusimos en el apartado 3 de las que realizaremos algunas simplificaciones: 𝜕(𝑈) 𝜕𝑡 +𝜕(𝐹) 𝜕𝑥 +𝜕(𝐺) 𝜕𝑦 +𝑆=0 ( 5-1) 𝑈=[𝜌 𝜌𝑢 𝜌𝑣 𝐸] ( 5-2) 𝐹= [ 𝜌𝑢 𝜌𝑢2+𝑝−𝜏𝑥𝑥 𝜌𝑢𝑣−𝜏𝑥𝑦 (𝐸+𝑝)𝑢−𝑢𝜏𝑥𝑥−𝑣𝜏𝑥𝑦+𝑞𝑥 ] ( 5-3) 80 Donde se aprecia con más claridad si cabe la dominancia de los mecanismos de transferencia de calor por convección en el caso de 𝑃𝑟=7 mientras que en el caso de 𝑃𝑟=0.7 ambos mecanismos son importantes. Si realizamos un análisis de órdenes de magnitud podemos estimar los espesores de la capa límite térmica. 𝛿𝑡 𝐷=1 √𝑅𝑒𝑃𝑟 ( 5-34) Teniendo en cuenta que ambos caos se han obtenido para un mismo número de Reynolds podemos obtener la siguiente relación: 𝛿𝑡0.7 𝛿𝑡7=√7 0.7=√10~3.16 ( 5-35) Esto se puede corroborar observando los perfiles de temperatura representados anteriormente: 𝛿𝑡0.7 𝛿𝑡7=2.5 0.9375=2.668 ( 5-36) Vemos como ambos números son del mismo orden de magnitud. No obstante, es conveniente comentar que en el caso de 𝑃𝑟=0.7 los efectos perturbadores de la pared del conducto aparecen ante de que se pueda desarrollar completamente el perfil de temperaturas. En la página siguiente mostramos el código MATLAB utilizado para la resolución de estos problemas. 81 %=======FLUJO ALREDEDOR DE UN CUADRADO CON TEMPERATURA=========% %LAS VARIABLES P_ SE CORRESPONDEN AL TERMINO PREDICTOR %==============================================================% close all;clear all; %========================PARÁMETROS============================% nx=319;%intervalos de partición en x ny=179;%intervalos de partición en y nt=35000;%intervalos de partición en t D=1; %ancho cavidad Dt=0.003; %incremento temporal Dx=35*D/nx; %incremento eje x Dy=6*D/ny; %incremento eje y %=======================NODOS CUADRADO==========================% i1=round(15/35*(nx+1)); i2=round(16/35*(nx+1)); j1=round(2.5/6*(ny+1)); j2=round(3.5/6*(ny+1)); %========================NÚMEROS ADIMENSIONALES=================% M=0.15; %número de Mach Re=100; %número de Reynolds gamma=1.4; dt=0.35; %incremento temperatura %=========================TEMPERATURAS Y VELOCIDAD==========================% T2=1+dt;%temperatura cuadrado T1=1; %temperatura fluido U=M; %velocidad cavidad %==========================COEFICIENTES==================================% a1=Dt/Dx; a2=Dt/Dy; a3=Dt/Fr^2; a4=-Dt/Pr/Re/(gamma-1)/Dx^2; a52=-Dt/Pr/Re/(gamma-1)/Dy^2; a5=4*Dt/(3*Re*Dx^2); a6=Dt/(Re*Dy^2); a7=Dt/(Re*Dx^2); a8=4*Dt/(3*Re*Dy^2); a9=Dt/(12*Re*Dx*Dy); a10=2*(a5+a6); a11=2*(a7+a8); %================MALLADO============================% x=linspace(0,35,nx+1);y=linspace(0,6,ny+1);[X,Y]=meshgrid(x,y) ; %mallado %==================CONDICIONES INICIALES=============% P_u=ones(nx+1,ny+1)*U;%velocidad x en los puntos del mallado P_v=zeros(nx+1,ny+1);%velocidad y en los puntos del mallado P_u(i1:i2,j1:j2)=0; %velocidad x en el cuadrado P_v(i1:i2,j1:j2)=0; %velocidad y en el cuadrado P_T=ones(nx+1,ny+1)*T1; %campo de temperatura en cada punto del mallado equilibrio P_rho=ones(nx+1,ny+1);% densidad en los puntos del mallado P_rhou=P_rho.*P_u;% densidad*u en los puntos del mallado P_rhov=P_rho.*P_v;%densidad*v en los puntos del mallado P_rhou2=P_rhou.*P_u;% densidad*u^2 en los puntos del mallado P_rhov2=P_rhov.*P_v;%densidad*v^2 en los puntos del mallado 82 P_rhouv=P_rhou.*P_v;%densidad*uv en los puntos del mallado P_E=(P_T/gamma/(gamma-1)-0.5*(P_v.^2+P_u.^2)).*P_rho;%energía total en los puntos del mallado P_Eu=P_E.*P_u;%energía*u total en los puntos del mallado P_Ev=P_E.*P_v;%energía*v total en los puntos del mallado u=ones(nx+1,ny+1)*U;%velocidad x en los puntos del mallado v=zeros(nx+1,ny+1);%velocidad y en los puntos del mallado u(i1:i2,j1:j2)=0; %velocidad x en el cuadrado v(i1:i2,j1:j2)=0;%velocidad v en el cuadrado T=ones(nx+1,ny+1)*T1; %campo de temperatura en cada punto del mallado equilibrio rho=ones(nx+1,ny+1);% densidad en los puntos del mallado rhou=rho.*u;% densidad*u en los puntos del mallado rhov=rho.*v;%densidad*v en los puntos del mallado rhou2=rhou.*u;% densidad*u^2 en los puntos del mallado rhov2=rhov.*v;%densidad*v^2 en los puntos del mallado rhouv=rhou.*v;%densidad*uv en los puntos del mallado E=(T/gamma/(gamma-1)-0.5*(v.^2+u.^2)).*rho;%energía total en los puntos del mallado Eu=E.*u;%energía*u total en los puntos del mallado Ev=E.*v;%energía*v total en los puntos del mallado %==================BUCLE TEMPORAL===========================% cont=1; k=1; for q=1:nt u(:,:)=rhou(:,:)./rho(:,:); v(:,:)=rhov(:,:)./rho(:,:); k=1+k; %===========================PREDICTOR=============================% pdgamma=(gamma-1)*(E-0.5*rho.*(u.^2+v.^2));%VARIABLE AUXILIAR PRESIÓN P_rho(2:nx,2:ny)=rho(2:nx,2:ny)-a1*(rhou(3:nx+1,2:ny)-rhou(2:nx,2:ny))- a2*(rhov(2:nx,3:ny+1)-rhov(2:nx,2:ny)); P_rhou(2:nx,2:ny)=rhou(2:nx,2:ny)-a1*(gamma-1)*(E(3:nx+1,2:ny)- E(2:nx,2:ny))+a1*0.5*(gamma-1)*((rhov2(3:nx+1,2:ny)+rhou2(3:nx+1,2:ny))- (rhov2(2:nx,2:ny)+rhou2(2:nx,2:ny)))-a1*(rhou2(3:nx+1,2:ny)- rhou2(2:nx,2:ny))-a2*(rhouv(2:nx,3:ny+1)-rhouv(2:nx,2:ny))- a10*u(2:nx,2:ny)+a5*(u(3:nx+1,2:ny)+u(1:nx1,2:ny))+a6*(u(2:nx,3:ny+1)+u(2:nx,1:ny-1))+a9*(v(3:nx+1,3:ny+1)+v(1:nx1,1:ny-1)-v(1:nx-1,3:ny+1)-v(3:nx+1,1:ny-1)); P_rhov(2:nx,2:ny)=rhov(2:nx,2:ny)-a2*(gamma-1)*(E(2:nx,3:ny+1)- E(2:nx,2:ny))+a1*0.5*(gamma-1)*((rhov2(2:nx,3:ny+1)+rhou2(2:nx,3:ny+1))- (rhov2(2:nx,2:ny)+rhou2(2:nx,2:ny)))-a1*((rhouv(3:nx+1,2:ny))- (rhouv(2:nx,2:ny)))-a2*(rhov2(2:nx,3:ny+1)-rhov2(2:nx,2:ny))- a11*v(2:nx,2:ny)+a7*(v(3:nx+1,2:ny)+v(1:nx1,2:ny))+a8*(v(2:nx,3:ny+1)+v(2:nx,1:ny-1))+a9*(u(3:nx+1,3:ny+1)+u(1:nx1,1:ny-1)-u(1:nx-1,3:ny+1)-u(3:nx+1,1:ny-1)); P_E(2:nx,2:ny)=E(2:nx,2:ny)- a1*(u(3:nx+1,2:ny).*(E(3:nx+1,2:ny)+pdgamma(3:nx+1,2:ny))- u(2:nx,2:ny).*(E(2:nx,2:ny)+pdgamma(2:nx,2:ny)))-a4*(T(3:nx+1,2:ny)- 2*T(2:nx,2:ny)+T(1:nx-1,2:ny))-a52*(T(2:nx,3:ny+1)- 2*T(2:nx,2:ny)+T(2:nx,1:ny-1))- a2*(v(2:nx,3:ny+1).*(E(2:nx,3:ny+1)+pdgamma(2:nx,3:ny+1))- v(2:nx,2:ny).*(E(2:nx,2:ny)+pdgamma(2:nx,2:ny))); 83 %=============================DECODIFICACIÓN VARIABLES==================% P_u=P_rhou./P_rho; P_v=P_rhov./P_rho; Vmod2(2:nx,2:ny)=P_v(2:nx,2:ny).^2+P_u(2:nx,2:ny).^2;%variable auxiliar módulo velocidad al cuadrado P_T(2:nx,2:ny)=(P_E(2:nx,2:ny)./P_rho(2:nx,2:ny)- 0.5*Vmod2(2:nx,2:ny))*gamma*(gamma-1); %========================CONDICIONES DE CONTORNO PREDICTOR==============% P_rho(1,1:ny+1)=rho(1,1:ny+1)-(a1/2)*(-rhou(3,1:ny+1)+4*rhou(2,1:ny+1)- 3*rhou(1,1:ny+1));%%%%x=0 densidad P_rho(nx+1,1:ny+1)=rho(nx+1,1:ny+1)+(a1/2)*(-rhou(nx1,1:ny+1)+4*rhou(nx,1:ny+1)-3*rhou(nx+1,1:ny+1));%%%x=d densidad P_rho(2:nx,1)=rho(2:nx,1)-(a2/2)*(-rhov(2:nx,3)+4*rhov(2:nx,2)- 3*rhov(2:nx,1))-(a1*U/2)*(rho(3:nx+1,1)-rho(1:nx-1,1));%%%%y=0 densidad P_rho(2:nx,ny+1)=rho(2:nx,ny+1)+a2/2*(-rhov(2:nx,ny-1)+4*rhov(2:nx,ny)- 3*rhov(2:nx,ny+1))-(a1*U/2)*(rho(3:nx+1,ny+1)-rho(1:nx-1,ny+1));%%%y=d densidad P_rho(i1,j1:j2)=rho(i1,j1:j2)+(a1/2)*(-rhou(i1-2,j1:j2)+4*rhou(i11,j1:j2)-3*rhou(i1,j1:j2));%%%%x=0 densidad cubo P_rho(i2,j1:j2)=rho(i2,j1:j2)-(a1/2)*(- rhou(i2+2,j1:j2)+4*rhou(i2+1,j1:j2)-3*rhou(i2,j1:j2));%%%x=d densidad cubo P_rho(i1:i2,j1)=rho(i1:i2,j1)+a2/2*(-rhov(i1:i2,j1-2)+4*rhov(i1:i2,j1-1)- 3*rhov(i1:i2,j1));%%%%y=0 densidad cubo P_rho(i1:i2,j2)=rho(i1:i2,j2)-(a2/2)*(- rhov(i1:i2,j2+2)+4*rhov(i1:i2,j2+1)-3*rhov(i1:i2,j2));%%%y=d densidad cubo P_u(2:nx,ny+1)=U;%%%%velocidad u y=d P_v(2:nx,ny+1)=0;%%%%velocidad v y=d P_u(2:nx,1)=U;%%%%velocidad u y=0 P_v(2:nx,1)=0;%%%%velocidad v y=0 P_u(nx+1,:)=U;%%%%velocidad u x=d P_v(nx+1,:)=0;%%%%velocidad v x=d P_u(1,:)=U;%%%%velocidad u x=0 P_v(1,:)=0;%%%%velocidad v x=0 P_u(i1:i2,j2)=0;%%%%velocidad u y=d cubo P_v(i1:i2,j2)=0;%%%%velocidad v y=d cubo P_u(i1:i2,j1)=0;%%%%velocidad u y=0 cubo P_v(i1:i2,j1)=0;%%%%velocidad v y=0 cubo P_u(i2,j1:j2)=0;%%%%velocidad u x=d cubo P_v(i2,j1:j2)=0;%%%%velocidad v x=d cubo P_u(i1,j1:j2)=0;%%%%velocidad u x=0 cubo P_v(i1,j1:j2)=0;%%%%velocidad v x=0 cubo P_u(i1:i2,j1:j2)=0; P_v(i1:i2,j1:j2)=0; P_T(:,1)=(4*P_T(:,2)-P_T(:,3))/3; P_T(:,ny+1)=(4*P_T(:,ny)-P_T(:,ny-1))/3; P_T(1,2:ny)=T1 ; P_T(nx+1,2:ny)=P_T(nx,2:ny)+(P_T(nx,2:ny)-P_T(nx-1,2:ny)); P_T(i1:i2,j1:j2)=T2; 84 P_E(:,1)=P_T(:,1).*P_rho(:,1)/gamma/(gamma-1); P_E(:,ny+1)=P_T(:,ny+1).*P_rho(:,ny+1)/gamma/(gamma-1); P_E(1,2:ny)=P_T(1,2:ny).*P_rho(1,2:ny)/gamma/(gamma-1); P_E(nx+1,2:ny)=P_T(nx+1,2:ny).*P_rho(nx+1,2:ny)/gamma/(gamma-1); P_rhou(2:nx,ny+1)=P_rho(2:nx,1+ny).* P_u(2:nx,ny+1);%%%%%%%densidad*u en y=d P_rhou(2:nx,1)=P_rho(2:nx,1).* P_u(2:nx,1);%%%%%%%densidad*u en y=0 P_rhou(1,:)=P_rho(1,:).* P_u(1,:);%%%%%%%densidad*u en x=0 P_rhou(nx+1,:)=P_rho(nx+1,:).* P_u(nx+1,:);%%%%%%%densidad*u en x=d P_rhov(2:nx,ny+1)=P_rho(2:nx,1+ny).* P_v(2:nx,ny+1);%%%%%%%densidad*v en y=d P_rhov(2:nx,1)=P_rho(2:nx,1).* P_v(2:nx,1);%%%%%%%densidad*v en y=0 P_rhov(1,:)=P_rho(1,:).* P_v(1,:);%%%%%%%densidad*v en x=0 P_rhov(nx+1,:)=P_rho(nx+1,:).* P_v(nx+1,:);%%%%%%%densidad*v en x=d %=======================ACTUALIZACIÓN DE VARIABLES=============% Vmod2(2:nx,2:ny)=P_v(2:nx,2:ny).^2+P_u(2:nx,2:ny).^2; P_T(2:nx,2:ny)=(P_E(2:nx,2:ny)./P_rho(2:nx,2:ny)- 0.5*Vmod2(2:nx,2:ny))*gamma*(gamma-1); P_rhou2=P_rhou.*P_u; P_rhouv=P_rhou.*P_v; P_rhov2=P_rhov.*P_v; P_Eu=P_E.*P_u; P_Ev=P_E.*P_v; %============================CORRECTOR========================% pdgamma=0.5*(pdgamma+(gamma-1).*(P_E-0.5.*P_rho.*(P_u.^2+P_v.^2))); rho(2:nx,2:ny)=0.5*(P_rho(2:nx,2:ny)+rho(2:nx,2:ny)-a1*(rhou(2:nx,2:ny)- rhou(1:nx-1,2:ny))-a2*(rhov(2:nx,2:ny)-rhov(2:nx,1:ny-1))); rhou(2:nx,2:ny)=0.5*(P_rhou(2:nx,2:ny)+rhou(2:nx,2:ny)-a1*(gamma1)*(P_E(2:nx,2:ny)-P_E(1:nx-1,2:ny))+a1*0.5*(gamma1)*(P_rhov2(2:nx,2:ny)+P_rhou2(2:nx,2:ny)-(P_rhov2(1:nx1,2:ny)+P_rhou2(1:nx-1,2:ny)))-a1*(rhou2(2:nx,2:ny)-rhou2(1:nx-1,2:ny))- a2*(P_rhouv(2:nx,2:ny)-P_rhouv(2:nx,1:ny-1))- a10*P_u(2:nx,2:ny)+a5*(P_u(3:nx+1,2:ny)+P_u(1:nx1,2:ny))+a6*(P_u(2:nx,3:ny+1)+P_u(2:nx,1:ny1))+a9*(P_v(3:nx+1,3:ny+1)+P_v(1:nx-1,1:ny-1)-P_v(1:nx-1,3:ny+1)- P_v(3:nx+1,1:ny-1))); rhov(2:nx,2:ny)=0.5*(P_rhov(2:nx,2:ny)+rhov(2:nx,2:ny)-a2*(gamma1)*(P_E(2:nx,2:ny)-P_E(2:nx,1:ny-1))+a2*0.5*(gamma1)*(P_rhov2(2:nx,2:ny)+P_rhou2(2:nx,2:ny)-(P_rhov2(2:nx,1:ny1)+P_rhou2(2:nx,1:ny-1)))-a1*((P_rhouv(2:nx,2:ny))-(P_rhouv(1:nx-1,2:ny)))- a2*(rhov2(2:nx,2:ny)-rhov2(2:nx,1:ny-1))- a11*v(2:nx,2:ny)+a7*(v(3:nx+1,2:ny)+v(1:nx1,2:ny))+a8*(v(2:nx,3:ny+1)+v(2:nx,1:ny-1))+a9*(u(3:nx+1,3:ny+1)+u(1:nx1,1:ny-1)-u(1:nx-1,3:ny+1)-u(3:nx+1,1:ny-1))); E(2:nx,2:ny)=0.5*(P_E(2:nx,2:ny)+E(2:nx,2:ny)- a1*(u(2:nx,2:ny).*(E(2:nx,2:ny)+pdgamma(2:nx,2:ny))-u(1:nx-1,2:ny).*(E(1:nx1,2:ny)+pdgamma(1:nx-1,2:ny)))-a4*(T(3:nx+1,2:ny)-2*T(2:nx,2:ny)+T(1:nx1,2:ny))-a52*(T(2:nx,3:ny+1)-2*T(2:nx,2:ny)+T(2:nx,1:ny-1))- a2*(v(2:nx,2:ny).*(E(2:nx,2:ny)+pdgamma(2:nx,2:ny))-v(2:nx,1:ny1).*(E(2:nx,1:ny-1)+pdgamma(2:nx,1:ny-1)))); %================ACTUALIZACIÓN DE VARIABLES=====================% u=rhou./rho; v=rhov./rho; Vmod2(2:nx,2:ny)=v(2:nx,2:ny).^2+u(2:nx,2:ny).^2; 85 T(2:nx,2:ny)=(E(2:nx,2:ny)./rho(2:nx,2:ny)- 0.5*Vmod2(2:nx,2:ny))*gamma*(gamma-1); %================C.C======================================% T(:,1)=(4*T(:,2)-T(:,3))/3;%%Y=0 pared caliente T(:,ny+1)=(4*T(:,ny)-T(:,ny-1))/3;%%%%y=d pared fria T(1,2:ny)=T1 ; % T(nx+1,2:ny)=T(nx,2:ny)+(T(nx,2:ny)-T(nx1,2:ny)); T(i1:i2,j2)=T2;%%%%velocidad u y=d cubo T(i1:i2,j1)=T2;%%%%velocidad u y=0 cubo T(i2,j1:j2)=T2;%%%%velocidad u x=d cubo T(i1,j1:j2)=T2;%%%%velocidad u x=0 cubo T(i1:i2,j1:j2)=T2; rho(1,2:ny)=(P_rho(1,2:ny)+rho(1,2:ny)-a1/2*(- P_rhou(3,2:ny)+4*P_rhou(2,2:ny)-3*P_rhou(1,2:ny)))*0.5;%%%%x=0 rho(nx+1,2:ny)=(P_rho(nx+1,2:ny)+rho(1,2:ny)+a1/2*(-P_rhou(nx1,2:ny)+4*P_rhou(nx,2:ny)-3*P_rhou(nx+1,2:ny)))*0.5;%%%x=d rho(2:nx,1)=(P_rho(2:nx,1)+rho(2:nx,1)-a2/2*(- P_rhov(2:nx,3)+4*P_rhov(2:nx,2)-3*P_rhov(2:nx,1))-a1*U/2*(P_rho(3:nx+1,1)- P_rho(1:nx-1,1)))*0.5;%%%%y=0 rho(2:nx,1+ny)=(P_rho(2:nx,1+ny)+rho(2:nx,1+ny)+a2/2*(-P_rhov(2:nx,ny1)+4*P_rhov(2:nx,ny)-3*P_rhov(2:nx,ny+1))-a1*U/2*(P_rho(3:nx+1,ny+1)- P_rho(1:nx-1,ny+1)))*0.5;%%%y=d rho(i1,j1:j2)=(rho(i1,j1:j2)+P_rho(i1,j1:j2)+(a1/2)*(-P_rhou(i12,j1:j2)+4*P_rhou(i1-1,j1:j2)-3*P_rhou(i1,j1:j2)))*0.5;%%%%x=0 densidad cubo rho(i2,j1:j2)=(rho(i2,j1:j2)+P_rho(i2,j1:j2)-(a1/2)*(- P_rhou(i2+2,j1:j2)+4*P_rhou(i2+1,j1:j2)-3*P_rhou(i2,j1:j2)))*0.5;%%%x=d densidad cubo rho(i1:i2,j1)=(rho(i1:i2,j1)+P_rho(i1:i2,j1)+a2/2*(-P_rhov(i1:i2,j12)+4*P_rhov(i1:i2,j1-1)-3*P_rhov(i1:i2,j1)))*0.5;%%%%y=0 densidad cubo rho(i1:i2,j2)=(rho(i1:i2,j2)+P_rho(i1:i2,j2)-(a2/2)*(- P_rhov(i1:i2,j2+2)+4*P_rhov(i1:i2,j2+1)-3*P_rhov(i1:i2,j2)))*0.5;%%%y=d densidad cubo u(2:nx,ny+1)=U;%%%%velocidad u y=d v(2:nx,ny+1)=0;%%%%velocidad v y=d u(2:nx,1)=U;%%%%velocidad u y=0 v(2:nx,1)=0;%%%%velocidad v y=0 u(nx+1,:)=U;%%%%velocidad u x=d v(nx+1,:)=0;%%%%velocidad v x=d u(1,:)=U;%%%%velocidad u x=0 v(1,:)=0;%%%%velocidad v x=0 u(i1:i2,j2)=0;%%%%velocidad u y=d cubo v(i1:i2,j2)=0;%%%%velocidad v y=d cubo u(i1:i2,j1)=0;%%%%velocidad u y=0 cubo v(i1:i2,j1)=0;%%%%velocidad v y=0 cubo u(i2,j1:j2)=0;%%%%velocidad u x=d cubo v(i2,j1:j2)=0;%%%%velocidad v x=d cubo u(i1,j1:j2)=0;%%%%velocidad u x=0 cubo v(i1,j1:j2)=0;%%%%velocidad v x=0 cubo u(i1:i2,j1:j2)=0; v(i1:i2,j1:j2)=0; 86 E(:,1)=T(:,1).*rho(:,1)/gamma/(gamma-1); E(:,ny+1)=T(:,ny+1).*rho(:,ny+1)/gamma/(gamma-1); E(1,2:ny)=T(1,2:ny).*rho(1,2:ny)/gamma/(gamma-1); E(nx+1,2:ny)=T(nx+1,2:ny).*rho(nx+1,2:ny)/gamma/(gamma-1); rhou(2:nx,ny+1)=rho(2:nx,1+ny).* u(2:nx,ny+1);%%%%%%%densidad*u en y=d rhou(2:nx,1)=rho(2:nx,1).*u(2:nx,1);%%%%%%%densidad*u en y=0 rhou(1,:)=rho(1,:).* u(1,:);%%%%%%%densidad*u en x=0 rhou(nx+1,:)=rho(nx+1,:).*u(nx+1,:);%%%%%%%densidad*u en x=d rhov(2:nx,ny+1)=rho(2:nx,1+ny).*v(2:nx,ny+1);%%%%%%%densidad*v en y=d rhov(2:nx,1)=rho(2:nx,1).*v(2:nx,1);%%%%%%%densidad*v en y=0 rhov(1,:)=rho(1,:).*v(1,:);%%%%%%%densidad*v en x=0 rhov(nx+1,:)=rho(nx+1,:).*v(nx+1,:);%%%%%%%densidad*v en x=d %===================ACTUALIZACIÓN DE VARIABLES=====================% Eu=E.*u; Ev=E.*v; rhou=rho.*u;% densidad*u en los puntos del mallado rhov=rho.*v;%densidad*v en los puntos del mallado rhou2=rhou.*u; rhouv=rhou.*v; rhov2=rhov.*v; q; %=============REPRESENTACIÓN GRÁFICA==========================% if k==50 disp(max(max(abs(v)))); % plot(H) % pause q % figure % contourf(T') % Tmod=sqrt(v'.^2+u'.^2); % figure % contourf(Tmod) % figure % quiver(linspace(0,2,nx+1),linspace(0,1,ny+1),u'./Tmod,v'./Tmod) k=1 pel(:,:,cont)=T; cont=cont+1; end H(q)=max(max(abs(u))); end quiver(linspace(0,2,nx+1),linspace(0,1,ny+1),u'./sqrt(v'.^2+u'.^2),v'./sqrt( v'.^2+u'.^2)) 87 6 APLICACIÓN: CAVIDAD RECTANGULAR CON DIFERENCIA DE TEMPERATURA EN LAS TAPAS El objetivo de esta aplicación será reproducir los patrones de convección natural originados por la inestabilidad que generan las fuerzas de flotabilidad en un fluido en reposo. Si tenemos un fluido en reposo y se introduce una pequeña perturbación que provoque un ligero desplazamiento vertical de la partícula, la densidad media de la partícula será menor que la de su entorno y tendrá una tendencia aún más ascendente. De forma opuesta si el desplazamiento de la partícula es negativo su densidad media será mayor que la de su entorno y esto provocará que su tendencia sea aún más descendente. Además de las fuerzas de flotabilidad aparecen otras fuerzas que tienden a amortiguar el efecto desestabilizador de la flotabilidad estas son los efectos de viscosidad y de conducción de calor. Muchos estudios de este problema han sido realizados mediante el uso de la aproximación de Boussinesq y en el caso incompresible, en este manual se trata el estudio con la resolución de las ecuaciones de Navier-Stokes completas utilizando un esquema de MacCormack, utilizado por K.V.Parchevsky [19] para realizar simulaciones a gran escala de la convección solar, por lo que se podría considerar como desarrollo futuro la aplicación de estas ecuaciones para la resolución de problemas de ese tipo. A continuación, vamos a realizar la formulación de la aplicación que aquí se va a resolver. Figura( 6-1 ) Se va a resolver el problema de convección de Rayleigh-Benard con la pared inferior a mayor temperatura, que como veremos es condición necesaria para que se genere la inestabilidad. Las paredes laterales vamos a considerarlas perfectamente adiabáticas. Comencemos con las ecuaciones que vamos a resolver. Partiremos de las ecuaciones de Navier-Stokes completas e iremos realizando algunas simplificaciones.  Se añadirán fuerzas másicas a las ecuaciones. 88 La velocidad característica se puede estimar de igualar en la ecuación de cantidad de movimiento el orden de magnitud de los términos convectivos y los de flotabilidad como: 𝑣𝑜=√Δ𝑇𝛽𝑔𝐻 ( 6-1) Además, el cuadrado del número de Froude es tal que: 𝐹𝑟2=𝑔𝐻Δ𝑇𝛽 𝑔𝐻 ~𝑂(1) ( 6-2) No podemos despreciar las fuerzas másicas lógicamente ya que son las encargadas de realizar el movimiento del fluido.  Los términos de disipación viscosa comparados con los de convección térmica se pueden expresar como: 𝜏:𝛻𝑣 𝜌𝑐𝑝𝑣 𝛻𝑇~𝜇𝑣02 𝐻2 𝜌𝑐𝑝𝑣0Δ𝑇 𝐻=𝜈 𝑐𝑝√𝑔𝐻Δ𝑇𝛽 𝐻Δ𝑇 ≪1 ( 6-3) Y por tanto podrán ser despreciados los efectos de disipación viscosa respecto a los de convección térmica. 𝜕(𝑈) 𝜕𝑡 +𝜕(𝐹) 𝜕𝑥 +𝜕(𝐺) 𝜕𝑦 +𝑆=0 ( 6-4) 𝑈=[𝜌 𝜌𝑢 𝜌𝑣 𝐸]=[𝑈1 𝑈2 𝑈3 𝑈4] ( 6-5) 𝐹=[𝜌𝑢 𝜌𝑢2+𝑝−𝜏𝑥𝑥 𝜌𝑢𝑣−𝜏𝑥𝑦 (𝐸+𝑝)𝑢+𝑞𝑥] ( 6-6) 𝐺= [ 𝜌𝑣 𝜌𝑢𝑣−𝜏𝑥𝑦 𝜌𝑣2+𝑝−𝜏𝑦𝑦 (𝐸+𝑝)𝑣+𝑞𝑦 ] ( 6-7) 𝑆=[00 𝜌𝑔 𝜌𝑣𝑔] ( 6-8) Para evitar valores muy altos de presión y temperatura y por tanto reducir el cálculo computacional, vamos a operar únicamente con las diferencias de temperatura y presión respecto a las variables en equilibrio hidrostático: 𝑝=𝑝𝑒𝑞+𝑝+ , 𝑇=𝑇𝑒𝑞+𝑇+, 𝐸=𝐸𝑒𝑞+𝐸+ ( 6-9) 89 Se introducen las siguientes variables adimensionales: 𝑥∗=𝑥𝐻, 𝑦∗=𝑦 𝐻,𝑣0=√ΔTβ𝑔𝐻=√ΔT∗𝑔𝐻 , 𝑢∗=𝑢 𝑣0 , 𝑣∗=𝑣 𝑣0, 𝑝∗ =𝑝+ 𝜌0v02, 𝑝𝑒𝑞 ∗=𝑝𝑒𝑞 𝜌0𝑅𝑔𝑇0, 𝜌∗=𝜌 𝜌0 , 𝑇∗=𝑇 𝑇0, 𝐸𝑒𝑞 ∗=𝐸𝑒𝑞 𝜌0𝑅𝑔𝑇0, 𝐸∗=𝐸+ 𝜌0v02 ( 6-10) Se definen los números adimenionales: 𝑅𝑒=𝜌0𝑣0𝐻 𝜇 , 𝑃𝑟=𝜇𝑐𝑝 𝑘=𝜈𝛼 , 𝐹𝑟2=ΔTβ=Δ𝑇 𝑇0=Δ𝑇∗,𝑀=𝑣0 𝑎0 ( 6-11) Que representan los Números de Reynolds, Prandtl, Froude y Mach respectivamente. Siendo 𝑣0 una velocidad característica, 𝜇 la viscosidad del fluido, 𝜌0 𝑦 𝑝0 la densidad y presión de referencia, 𝜅 la conductividad térmica del fluido, 𝑔 la aceleración gravitatoria, H una dimensión característica del problema y 𝑎0 la velocidad del sonido. Con estas nuevas variables adimensionales, de las que de ahora en adelante se suprimirá el asterisco, el sistema de ecuaciones resultante será: 𝑑(𝑈) 𝑑𝑡 +𝑑(𝐹) 𝑑𝑥 +𝑑(𝐺) 𝑑𝑦 +𝑆=0 ( 6-12) 𝑈=[𝜌 𝜌𝑢 𝜌𝑣 𝐸∗]=[𝑈1 𝑈2 𝑈3 𝑈4] ( 6-13) 𝐹=[𝜌𝑢 𝜌𝑢2+p∗−𝜏𝑥𝑥 𝜌𝑢𝑣−𝜏𝑥𝑦 (𝐸∗+p∗)𝑢+𝑞𝑥] ( 6-14) 𝐺= [ 𝜌𝑣 𝜌𝑢𝑣−𝜏𝑥𝑦 𝜌𝑣2+𝑝∗−𝜏𝑦𝑦 (𝐸∗+p∗)𝑣+𝑞𝑦 ] ( 6-15) 𝑆= [ 00 𝜌−𝜌𝑒𝑞 𝐹𝑟2 𝜌𝑣−𝛾 𝛾−1 𝜌𝑒𝑞𝑣 𝐹𝑟2+𝑝𝑒𝑞 (𝛾−1)𝑀2∇.𝑽 ] ( 6-16) 96 𝑅𝑎=659.58 Figura( 6-3 ) Figura( 6-4 ) 97 𝑅𝑎=4717 Figura( 6-5 ) Figura( 6-6 ) 98 Observamos como en el primer caso, con el número de Rayleigh inferior al crítico, el campo de temperaturas se muestra sin alterar comparado con la condición inicial. En el segundo caso el número de Rayleigh ha superado el umbral a partir del cual los efectos de flotabilidad superan a los viscosos y a los de conducción de calor y el fluido comienza a ascender por el centro y a descender por las paredes formando este patrón celular tan representativo. Este patrón de flujo provoca que el fluido caliente de la placa de abajo asciende por el centro y el fluido mas frio desciende por las paredes, provocando la distribución de temperaturas mostrada en la figura. Otro aspecto del problema de Rayleigh-Bénard que se puede estudiar es modificar las condiciones de contorno en las paredes. Figura( 6-7 ) Las paredes laterales vamos a considerarlas perfectamente conductoras. La situación es equivalente a una cavidad embebida en una plancha metálica con diferencia de temperatura entre caras. En la figura 6-7 se muestra el esquema general del problema. El problema se resuelva idénticamente al anterior a excepción de la condición de contorno en temperatura que en este caso se reduce a imponer que la temperatura asociada al movimiento sea nula em las paredes, de esta forma la temperatura en las paredes será la de equilibrio que será la distribución que tendrá la plancha metálica. 99 Figura( 6-8 ) Observamos como el efecto de las paredes conductoras tiene un efecto amortiguador para la convección ya que acentúa más los efectos de conducción de calor haciendo más difícil la aparición de un régimen estacionario. Cuando este aparece lo hace con este patrón diferente al anterior debido a que la distribución de temperaturas en las paredes tiene que ser lineal por la condición de contorno de paredes conductoras que se ha impuesto. Un estudio exhaustivo para determinar el número de Rayleigh crítico nos mostraría que para una relación de aspecto dada el número de Rayleigh crítico para el caso de paredes conductores sería mayor que para paredes adiabáticas ya que como hemos comentado esta condición de contorno acentúa los efectos de conducción de calor que amortiguan la inestabilidad producida por las fuerzas de flotabilidad. 100 %==============cavidad con diferencia de temperaturas==============% close all;clear all; %=======================Parámetros==========================% nx=50;%intervalos de partición en x ny=50;%intervalos de partición en y nt=150000 ;%intervalos de partición en t D=1; %ancho cavidad Dt=0.0001 ; %incremento temporal Dx=2*D/nx; %incremento eje x Dy=D/ny; %incremento eje y %=============NUMEROS ADIMENSIONALES==========================% dt=0.5; %INCREMENTO ADIMENSIONAL DE TEMPERATURA M=0.1; %NÚMERO DE MACH Re=30; %NÚMERO DE REYNOLDS Pr=0.73; %NÚMERO DE PRANDTL gamma=1.4; %COEFICIENTE DE EXPANSIÓN ADIABÁTICO Fr=sqrt(dt);%NÚMERO DE FROUDE %============================TEMPERATURAS===============================t T2=1-dt/2; T1=1+dt/2; %============================COEFICIENTES===============================% a1=Dt/Dx; a2=Dt/Dy; a3=Dt/Fr^2; a5=4*Dt/(3*Re*Dx^2); a6=Dt/(Re*Dy^2); a7=Dt/(Re*Dx^2); a8=4*Dt/(3*Re*Dy^2); a9=Dt/(12*Re*Dx*Dy); a10=2*(a5+a6); a11=2*(a7+a8); %===============MALLADO===================================% x=linspace(0,2,nx+1);y=linspace(0,1,ny+1);[X,Y]=meshgrid(x,y) ; %mallado %======================CONDICIONES INICIALES===============% P_u=10^-10*ones(nx+1,ny+1).*sin(2*pi*X).*cos(pi*Y);%velocidad residual direccion x en los puntos del mallado P_v=-10^-10*ones(nx+1,ny+1).*cos(2*pi*X).*sin(pi*Y);%velocidad residual dirección y en los puntos del mallado P_Teq=(linspace(T1,T2,ny+1)'*ones(1,nx+1))'; %campo de temperatura en cada punto del mallado equilibrio P_T=zeros(nx+1,ny+1); %campo de temperatura DEL MOVIMIENTO en cada punto del mallado equilibrio P_rhoeq=P_Teq.^(gamma*M^2/dt^2-1);% densidad en los puntos del mallado P_rho=P_rhoeq;% densidad en los puntos del mallado P_rhou=P_rho.*P_u;% densidad*u en los puntos del mallado P_rhov=P_rho.*P_v;%densidad*v en los puntos del mallado P_rhou2=P_rhou.*P_u;% densidad*u^2 en los puntos del mallado P_rhov2=P_rhov.*P_v;%densidad*v^2 en los puntos del mallado P_rhouv=P_rhou.*P_v;%densidad*uv en los puntos del mallado P_E=(P_T.*P_rho-(P_rhoeq-P_rho).*P_Teq)/gamma/(gamma1)/M^2+P_rhou2/2+P_rhov2/2; 101 P_Eu=P_E.*P_u; P_Ev=P_E.*P_v; u=10^-10*ones(nx+1,ny+1).*sin(2*pi*X).*cos(pi*Y);%velocidad residual direccion x en los puntos del mallado v=-10^-10*ones(nx+1,ny+1).*cos(2*pi*X).*sin(pi*Y);%velocidad residual dirección y en los puntos del mallado Teq=(linspace(T1,T2,ny+1)'*ones(1,nx+1))'; %campo de temperatura en cada punto del mallado equilibrio T=zeros(nx+1,ny+1); rhoeq=Teq.^(gamma*M^2/dt^2-1);% densidad en los puntos del mallado rho=rhoeq; rhou=rho.*u;% densidad*u en los puntos del mallado rhov=rho.*v;%densidad*v en los puntos del mallado rhou2=rhou.*u;% densidad*u^2 en los puntos del mallado rhov2=rhov.*v;%densidad*v^2 en los puntos del mallado rhouv=rhou.*v;%densidad*uv en los puntos del mallado E=(T.*rho-(rhoeq-rho).*Teq)/gamma/(gamma-1)/M^2+rhou2/2+rhov2/2; Eu=E.*u; Ev=E.*v; P_eq=rhoeq.*Teq; %PRESIÓN EQUILIBRIO Ra1=Re^2*Pr; %NÚMERO RAYLEIGH STANDAR Ra2=Re^2/Fr^2*Pr*(max(max(rho))-min(min(rho))); %NÚMERO DE RAYLEIGH CORREGIDO pause %========================BUCLE TEMPORAL===============================% k=0; cont=1; for q=1:nt k=k+1; %===========================================PREDICTOR======================= % P_rho(2:nx,2:ny)=rho(2:nx,2:ny)-a1*(rhou(3:nx+1,2:ny)-rhou(2:nx,2:ny))- a2*(rhov(2:nx,3:ny+1)-rhov(2:nx,2:ny)); P_rhou(2:nx,2:ny)=rhou(2:nx,2:ny)-a1*(gamma-1)*(E(3:nx+1,2:ny)- E(2:nx,2:ny)-rhou2(3:nx+1,2:ny)/2+rhou2(2:nx,2:ny)/2rhov2(3:nx+1,2:ny)/2+rhov2(2:nx,2:ny)/2)-a1*(rhou2(3:nx+1,2:ny)- rhou2(2:nx,2:ny))-a2*(rhouv(2:nx,3:ny+1)-rhouv(2:nx,2:ny))- a10*u(2:nx,2:ny)+a5*(u(3:nx+1,2:ny)+u(1:nx1,2:ny))+a6*(u(2:nx,3:ny+1)+u(2:nx,1:ny-1))+a9*(v(3:nx+1,3:ny+1)+v(1:nx1,1:ny-1)-v(1:nx-1,3:ny+1)-v(3:nx+1,1:ny-1)); P_rhov(2:nx,2:ny)=rhov(2:nx,2:ny)-a2*(gamma-1)*(E(2:nx,3:ny+1)- E(2:nx,2:ny)-rhou2(2:nx,3:ny+1)/2+rhou2(2:nx,2:ny)/2rhov2(2:nx,3:ny+1)/2+rhov2(2:nx,2:ny)/2)-a2*(rhov2(2:nx,3:ny+1)- rhov2(2:nx,2:ny))-a1*(rhouv(3:nx+1,2:ny)-rhouv(2:nx,2:ny))- a11*v(2:nx,2:ny)+a7*(v(3:nx+1,2:ny)+v(1:nx1,2:ny))+a8*(v(2:nx,3:ny+1)+v(2:nx,1:ny-1))+a9*(u(3:nx+1,3:ny+1)+u(1:nx1,1:ny-1)-u(1:nx-1,3:ny+1)-u(3:nx+1,1:ny-1))-(rho(2:nx,2:ny)- rhoeq(2:nx,2:ny))*Dt/(Fr^2); P_E(2:nx,2:ny)=E(2:nx,2:ny)-a1*(gamma)*(Eu(3:nx+1,2:ny)- Eu(2:nx,2:ny))+a1*(gamma1)*0.5*(u(3:nx+1,2:ny).*(rhou2(3:nx+1,2:ny)+rhov2(3:nx+1,2:ny))- u(2:nx,2:ny).*(rhou2(2:nx,2:ny)+rhov2(2:nx,2:ny)))+a1/(Dx*Pr*Re*M^2*(gamma1))*(T(3:nx+1,2:ny)-2*T(2:nx,2:ny)+T(1:nx-1,2:ny))- a2*(gamma)*(Ev(2:nx,3:ny+1)-Ev(2:nx,2:ny))+a2*0.5*(gamma- 102 1)*(v(2:nx,3:ny+1).*(rhou2(2:nx,3:ny+1)+rhov2(2:nx,3:ny+1))- v(2:nx,2:ny).*(rhou2(2:nx,2:ny)+rhov2(2:nx,2:ny)))+a2/(Dy*Pr*Re*M^2*(gamma1))*(T(2:nx,3:ny+1)-2*T(2:nx,2:ny)+T(2:nx,1:ny-1))-(rho(2:nx,2:ny)- rhoeq(2:nx,2:ny)*gamma/(gamma-1)).*v(2:nx,2:ny)*Dt/(Fr^2)-1/M^2/(gamma1)*P_eq(2:nx,2:ny).*(a1*(u(3:nx+1,2:ny)-u(2:nx,2:ny))+a2*(v(2:nx,3:ny+1)- v(2:nx,2:ny))); %===========================DECODIFICACIÓN=======================% P_u=P_rhou./P_rho; P_v=P_rhov./P_rho; Vmod2=P_v.^2+P_u.^2; P_T=(P_rhoeq-P_rho).*P_Teq./P_rho+(P_E./P_rho-Vmod2)*M^2*gamma*(gamma-1); %===========================C.C===================================% P_T(:,1)=0;%%Y=0 pared caliente P_T(:,ny+1)=0;%%%%y=d pared fria P_T(1,2:ny)=(4*P_T(2,2:ny)-P_T(3,2:ny))/3;%%%%%X=0 dT/dx=0 P_T(nx+1,2:ny)=(4*P_T(nx,2:ny)-P_T(nx-1,2:ny))/3;%%%%%X=d dT/dx=0 P_rho(1,:)=rho(1,1:ny+1)-(a1/2)*(-rhou(3,1:ny+1)+4*rhou(2,1:ny+1)- 3*rhou(1,1:ny+1));%%%%x=0 densidad P_rho(nx+1,:)=rho(nx+1,1:ny+1)+(a1/2)*(-rhou(nx1,1:ny+1)+4*rhou(nx,1:ny+1)-3*rhou(nx+1,1:ny+1));%%%x=d densidad P_rho(2:nx,1)=rho(2:nx,1)-(a2/2)*(-rhov(2:nx,3)+4*rhov(2:nx,2)- 3*rhov(2:nx,1));%%%%y=0 densidad P_rho(2:nx,ny+1)=rho(2:nx,ny+1)+a2/2*(-rhov(2:nx,ny-1)+4*rhov(2:nx,ny)- 3*rhov(2:nx,ny+1));%%%y=d densidad P_u(:,ny+1)=0;%%%%velocidad u y=d P_v(:,ny+1)=0;%%%%velocidad v y=d P_u(:,1)=0;%%%%velocidad u y=0 P_v(:,1)=0;%%%%velocidad v y=0 P_v(1,:)=0; P_v(nx+1,:)=0; P_u(1,:)=0;%%%%%X=0 du/dx=0 P_u(nx+1,:)=0;%%%%%X=d du/dx=0 %====================================ACTUALIZACIÓN DE VARIABLES===================% P_rhou=P_rho.*P_u;% densidad*u en los puntos del mallado P_rhov=P_rho.*P_v;%densidad*v en los puntos del mallado P_rhou2=P_rhou.*P_u;% densidad*u^2 en los puntos del mallado P_rhov2=P_rhov.*P_v;%densidad*v^2 en los puntos del mallado P_rhouv=P_rhou.*P_v;%densidad*uv en los puntos del mallado P_E=(P_T.*P_rho-(P_rhoeq-P_rho).*P_Teq)/gamma/(gamma1)/M^2+P_rhou2/2+P_rhov2/2; P_Eu=P_E.*P_u; P_Ev=P_E.*P_v; %================================CORRECTOR================================% 103 rho(2:nx,2:ny)=0.5*(P_rho(2:nx,2:ny)+rho(2:nx,2:ny)- a1*(P_rhou(2:nx,2:ny)-P_rhou(1:nx-1,2:ny))-a2*(P_rhov(2:nx,3:ny+1)- P_rhov(2:nx,1:ny-1))); rhou(2:nx,2:ny)=0.5*(rhou(2:nx,2:ny)+P_rhou(2:nx,2:ny)-a1*(gamma1)*(E(2:nx,2:ny)-E(1:nx-1,2:ny)-rhou2(2:nx,2:ny)/2+rhou2(1:nx-1,2:ny)/2rhov2(2:nx,2:ny)/2+rhov2(1:nx-1,2:ny)/2)-a1*(P_rhou2(2:nx,2:ny)- P_rhou2(1:nx-1,2:ny))-a2*(P_rhouv(2:nx,2:ny)-P_rhouv(2:nx,1:ny-1))- a10*P_u(2:nx,2:ny)+a5*(P_u(3:nx+1,2:ny)+P_u(1:nx1,2:ny))+a6*(P_u(2:nx,3:ny+1)+P_u(2:nx,1:ny1))+a9*(P_v(3:nx+1,3:ny+1)+P_v(1:nx-1,1:ny-1)-P_v(1:nx-1,3:ny+1)- P_v(3:nx+1,1:ny-1))); rhov(2:nx,2:ny)=0.5*(rhov(2:nx,2:ny)+P_rhov(2:nx,2:ny)-a2*(gamma1)*(E(2:nx,2:ny)-E(2:nx,1:ny-1)-rhou2(2:nx,2:ny)/2+rhou2(2:nx,1:ny-1)/2rhov2(2:nx,2:ny)/2+rhov2(2:nx,1:ny-1)/2)-a2*(P_rhov2(2:nx,2:ny)- P_rhov2(2:nx,1:ny-1))-a1*(P_rhouv(2:nx,2:ny)-P_rhouv(1:nx-1,2:ny))- a11*P_v(2:nx,2:ny)+a7*(P_v(3:nx+1,2:ny)+P_v(1:nx1,2:ny))+a8*(P_v(2:nx,3:ny+1)+P_v(2:nx,1:ny1))+a9*(P_u(3:nx+1,3:ny+1)+P_u(1:nx-1,1:ny-1)-P_u(1:nx-1,3:ny+1)- P_u(3:nx+1,1:ny-1))-(P_rho(2:nx,2:ny)-P_rhoeq(2:nx,2:ny))*Dt/(Fr^2)); E(2:nx,2:ny)=0.5*(E(2:nx,2:ny)+P_E(2:nx,2:ny)- a1*(gamma)*(P_Eu(2:nx,2:ny)-P_Eu(1:nx-1,2:ny))+a1*(gamma1)*0.5*(P_u(2:nx,2:ny).*(P_rhou2(2:nx,2:ny)+P_rhov2(2:nx,2:ny))-P_u(1:nx1,2:ny).*(P_rhou2(1:nx-1,2:ny)+P_rhov2(1:nx1,2:ny)))+a1/(Dx*Pr*Re*M^2*(gamma-1))*(P_T(3:nx+1,2:ny)- 2*P_T(2:nx,2:ny)+P_T(1:nx-1,2:ny))-a2*(gamma)*(P_Ev(2:nx,2:ny)- P_Ev(2:nx,1:ny-1))+a2*(gamma1)*0.5*(P_v(2:nx,2:ny).*(P_rhou2(2:nx,2:ny)+P_rhov2(2:nx,2:ny))- P_v(2:nx,1:ny-1).*(P_rhou2(2:nx,1:ny-1)+P_rhov2(2:nx,1:ny1)))+a2/(Dy*Pr*Re*M^2*(gamma-1))*(P_T(2:nx,3:ny+1)- 2*P_T(2:nx,2:ny)+P_T(2:nx,1:ny-1))-(P_rho(2:nx,2:ny)- P_rhoeq(2:nx,2:ny)*gamma/(gamma-1)).*P_v(2:nx,2:ny)*Dt/(Fr^2)-1/M^2/(gamma1)*P_eq(2:nx,2:ny).*(a1*(P_u(2:nx,2:ny)-P_u(1:nx1,2:ny))+a2*(P_v(2:nx,2:ny)-P_v(2:nx,1:ny-1)))); %===========================DECODIFICACIÓN=======================% u=rhou./rho; v=rhov./rho; Vmod2=v.^2+u.^2; T=(rhoeq-rho).*Teq./rho+(E./rho-Vmod2)*M^2*gamma*(gamma-1); %===========================C.C===================================% T(:,1)=0;%%Y=0 pared caliente T(:,ny+1)=0;%%%%y=d pared fria T(1,2:ny)=(4*T(2,2:ny)-T(3,2:ny))/3;%%%%%X=0 dT/dx=0 T(nx+1,2:ny)=(4*T(nx,2:ny)-T(nx-1,2:ny))/3;%%%%%X=d dT/dx=0 rho(1,:)=(P_rho(1,:)+rho(1,:)-a1/2*(-P_rhou(3,:)+4*P_rhou(2,:)- 3*P_rhou(1,:)))*0.5;%%%%x=0 rho(nx+1,:)=(P_rho(nx+1,:)+rho(1,:)+a1/2*(-P_rhou(nx1,:)+4*P_rhou(nx,:)-3*P_rhou(nx+1,:)))*0.5;%%%x=d rho(2:nx,1)=(P_rho(2:nx,1)+rho(2:nx,1)-a2/2*(- P_rhov(2:nx,3)+4*P_rhov(2:nx,2)-3*P_rhov(2:nx,1)))*0.5;%%%%y=0 rho(2:nx,1+ny)=(P_rho(2:nx,1+ny)+rho(2:nx,1+ny)+a2/2*(-P_rhov(2:nx,ny1)+4*P_rhov(2:nx,ny)-3*P_rhov(2:nx,ny+1)))*0.5;%%%y=d u(:,ny+1)=0;%%%%velocidad u y=d v(:,ny+1)=0;%%%%velocidad v y=d u(:,1)=0;%%%%velocidad u y=0 104 v(:,1)=0;%%%%velocidad v y=0 v(1,:)=0;%%%%%X=0 dv/dx=0 v(nx+1,:)=0;%%%%%X=d dv/dx=0 u(1,:)=0;%%%%%X=0 du/dx=0 u(nx+1,:)=0;%%%%%X=d du/dx=0 %====================================ACTUALIZACIÓN DE VARIABLES===================% rhou=rho.*u;% densidad*u en los puntos del mallado rhov=rho.*v;%densidad*v en los puntos del mallado rhou2=rhou.*u;% densidad*u^2 en los puntos del mallado rhov2=rhov.*v;%densidad*v^2 en los puntos del mallado rhouv=rhou.*v;%densidad*uv en los puntos del mallado E=(T.*rho-(rhoeq-rho).*Teq)/gamma/(gamma-1)/M^2+rhou2/2+rhov2/2; Eu=E.*u; Ev=E.*v; %==================REPRESENTACIÓN GRÁFICA Y CONVERGENCIA===================% if k==500 Dt; q; pel(:,:,cont)=T; cont=cont+1; k=0; end H(q)=max(max(abs(v))); if max(max(u))>=20 q=nt; pause end end quiver(linspace(0,2,nx+1),linspace(0,1,ny+1),u'./sqrt(v'.^2+u'.^2),v'./sqrt( v'.^2+u'.^2)) 105 7 CONCLUSIONES Y DESARROLLOS FUTUROS El objetivo principal de este proyecto ha sido demostrar que, hoy en día, los estudiantes de disciplinas como la Dinámica de Fluidos y la Transferencia de Calor pueden resolver eficientemente problemas realistas haciendo uso de conocimientos básicos de programación y de cálculo numérico adquiridos en asignaturas cursadas en años anteriores. El uso del ordenador permite al alumno un mejor entendimiento de las leyes fundamentales de la fluidodinámica, sin que tenga que recurrir a manipulaciones engorrosas ni a simplificaciones que hagan perder al problema todo su interés, buscando una solución las ecuaciones que gobiernan el problema. Para ilustrar esto se ha llevado a cabo un análisis numérico mediante el método de MacCormack aplicado a diversos problemas fluidodinámicos como el flujo alrededor de un cilindro cuadrado que viaja a velocidad constante por un conducto, la capa límite térmica alrededor del mismo y el régimen laminar no lineal de la convección natural de RayleighBénard en el interior de una cavidad rectangular. En el capítulo 3 se ha expuesto el método de MacCormack explicito, el esquema multi-etapa de diferencias finitas, de segundo orden tanto el tiempo como en el espacio, que hemos utilizado para resolver los sistemas de ecuaciones que se han presentado en las diferentes aplicaciones. Se han desarrollado dos ejemplos, el primero unidimensional, la resolución de la ecuación lineal de onda, y el segundo, ya más enfocado a la resolución del trabajo, un caso bidimensional en el cual la convección esta llevada por la tapadera, utilizaremos estos ejemplos como validadores del método ya que como podemos observar el acuerdo con los resultados de la bibliografía es excepcional. En el capítulo 4, se ha estudiado el flujo compresible alrededor de un cilindro cuadrado, primero realizamos una comparación de los coeficientes de sustentación y resistencia con los obtenidos en la literatura, obteniendo un acuerdo excelente. Después analizamos como afectaba la relación de aspecto del dominio computacional a los resultados para buscar una solución de compromiso entre gasto computacional y resultados que se aproximasen a un cilindro cuadrado inmerso en una corriente incidente. Una vez seleccionado el dominio computacional estudiamos el fenómeno de ‘vortex shedding’ y la frecuencia con la que los torbellinos se emitían en función del número de Reynolds, el acuerdo con los datos de la bibliografía es muy bueno teniendo en cuenta la limitación computacional que tenemos. El incremento temporal que nos exige la condición CFL es muy restrictivo y el espaciamiento de la malla se vuelve restrictivo cuando aumentamos el número de Reynolds lo que penaliza aún más la condición CFL. Dos posibles soluciones para esto podrían ser un refinamiento de malla selectivo en zonas donde sea necesario, en este caso en el entorno del cuadrado se podría hacer una malla más fina, y en el entorno de las secciones de entrada y salida y paredes laterales se podría tener una malla más gruesa ya que lo único que queremos en esos puntos es imponer condiciones de contorno y el desarrollo del método de MacCormack modificado que se explica brevemente en el capítulo tres que relaja la condición CFL al dividir el proceso en etapas unidimensionales. En el capítulo 5, se ha desarrollado una breve aplicación como primer contacto a la inclusión de la temperatura en las ecuaciones. Los perfiles de temperaturas obtenidos eran consecuentes con lo esperado y la estimación de órdenes de magnitud del espesor de la capa límite térmica también lo eran. Este problema requiere una capacidad computacional elevada por eso no se