88 Equa ion Chap e 1 Sec ion 1
T abajo Fin de G ado
G ado en Ingenie ía Ae oespacial
Aplicaciones del mé odo de MacCo mack a di e sos
p oblemas luidomecánicos
Au o : Manuel Ca eño Ruiz
Tu o : Miguel Pé ez-Sabo id Sánchez-Pas o
Dep. Ingenie ía Ae oespacial y Mecánica de Fluidos
Escuela Técnica Supe io de Ingenie ía
Uni e sidad de Se illa
Se illa, 2016
2
3
T abajo Fin de G ado
G ado en Ingenie ía Ae oespacial
Aplicaciones del mé odo de MacCo mack a
di e sos p oblemas luidomecánicos
Au o :
Manuel Ca eño Ruiz
Tu o :
Miguel Pé ez-Sabo id Sánchez-Pas o
P o eso i ula
Dep.Ingenie ía Ae oespacial y Mecánica de Fluidos
Escuela Técnica Supe io de Ingenie ía
Uni e sidad de Se illa
Se illa, 2016
4
5
T abajo Fin de ado: Aplicaciones del mé odo de MacCo mack a di e sos p oblemas luidomecánicos
Au o :
Manuel Ca eño Ruiz
Tu o :
Miguel Pé ez-Sabo id Sánchez-Pas o
El ibunal nomb ado pa a juzga el P oyec o a iba indicado, compues o po los siguien es miemb os:
P esiden e:
Vocales:
Sec e a io:
Acue dan o o ga le la cali icación de:
Se illa, 2016

6
El Sec e a io del T ibunal
7
A mi amilia,
amigos y compañe os.
8
Ag adecimien os
Deseo exp esa mi más p o undo ag adecimien o a odos los docen es con los que he enido el place de
eco e es e du o camino que aho a cie o con es e documen o, especialmen e a Don Miguel Pé ez-Sabo id
Sánchez-Pas o po el ines imable apoyo ecibido du an e es os meses de abajo.
También me gus a ía ag adece a mi amilia y amigos el apoyo que me han b indado du an e es os cua o
años, sin ningún asomo de duda sin ellos hoy no es a ía aquí.
Po úl imo, me gus a ía ag adece a mis compañe os la ayuda ecibida, especialmen e a Ma ía Paz Luque
Alcán a a y Ca men Vela de López de Ayala, cuya compañía du an e la gos días de abajo han
hecho que odo sea mucho más sencillo y han con ibuido de mane a inigualable a mi o mación an o
pe sonal como uni e si a ia.
Manuel Ca eño Ruiz
Se illa, 2016
9
Resumen
Es as memo ias ecogen unas ins ucciones y pasos sencillos median e los cuales el lec o pueda se capaz de
ene una comp ensión del mé odo de MacCo mack que le pe mi i á implemen a lo en un so wa e de
cálculo numé ico como MATLAB o simila pa a esol e di e sos p oblemas luidodinámicos.
Además de anliza la alidez del mé odo compa ándo con esul ados de la li e a u a se p ocede á a analiza
los esul ados de p oblemas como el lujo al ededo de un cilind o cuad ado, la con ección en una ca idad
cuad ada lle ada po la apade a o el p oblema de con ección na u al de Rayleigh-Béna d.
Los dos p oblemas es udiados en más p o undidad son la emisión de o bellinos que se p oduce a pa i un
núme o de Reynolds cuando una co ien e ebo dea un obs áculo, en es e caso un cuad ado, y algunos
aspec os elacionados con el p oblema de con ección de Rayleigh-Bena d, como el núme o de Reynolds
c í ico o la in luencia de las condiciones de con o no en las pa edes de la ca idad.
Como aspec o di e enciado a la mayo ía de esul ados de la bibliog a ía pa a el p oblema de Rayleigh-
Béna d en es a memo ia se con empla la esolución de las ecuaciones de Na ie -S okes comp esibles, sin
aplica la ap oximación de Boussinesq lo que nos pe mi á comp oba la alidez de es a.
16
Hen i Béna d lle ó a cabo una se ie de expe imen os en capas de luidos de ele ada conduc i idad é mica
calen adas desde abajo y cuya su supe icie supe io es aba expues a al ai e del ambien e, lo que dio
comienzo al es udio de los enómenos de con ección na u al.
Bé na d pudo obse a la exis encia de un mo imien o en el luido que alcanzaba un égimen es aciona io,
as un égimen ansi o io inicial, donde se dis inguía en la supe icie una es uc u a similia a la de un panal
de abejas de celdas hexagonales egula es como se obse a en la igu a 1-4. Realizó di e sas in es igaciones
del lujo en el in e io de la celda, las cuales e ela on que el luido ascendía po el cen o de la celda y
descendía a lo la go del pe íme o hexagonal.
Figu a( 1-4 )
El p ime in es igado que desa olló una eo ía que ponía de mani ies o de o ma cla a los mecanismos
ísicos in oluc ados en el enómeno de la con ección na u al, ue Lo d Rayleigh.
Rayleigh, esal ó que además de las ue zas de lo abilidad, deben ene se en cuen a, cuando un luido cuya
densidad no es uni o me se encuen a en p esencia de la g a edad; los mecanismos que ienden a
con a es a su e ec o: la icción, o iginada po las ue zas de iscosidad, y el de conducción de calo , que
iende a homogeneiza en el seno del luido el campo de empe a u as y, po an o, ambién el de densidades.
Rayleigh conside ó una capa de luido en e dos placas pa alelas in ini as a di e en e empe a u a dispues as
pe pendicula men e a la di ección de la g a edad, con la placa in e io a mayo empe a u a. Rayleigh
linealizo las ecuaciones de Na ie -S okes en o no a la si uación de equilib io, descub iendo así que en el
sis ema lineal de ecuaciones en de i adas pa ciales apa ece el pa áme o adimensional, denominado núme o
de Rayleigh:
𝑅𝑎=𝑔𝛽Δ𝑇𝐻3
𝜈𝛼
( 1-1)
donde β es el coe icien e de dila ación é mica, ν la iscosidad cinemá ica, α la di usi idad é mica del luido,
g es la acele ación de la g a i a o ia, H es una longi ud ca ac e ís ica del p oblema y ΔT es la di e encia de
empe a u a en e ambas placas.
El núme o de Rayleigh es el Núme o adimensional que mide la impo ancia ela i a en e los e ec os de las
ue zas de lo abilidad y los e ec os iscosos y de conducción é mica. Exis e un cie o alo umb al del
núme o de Rayleigh a pa i del cual los e ec os de lo abilidad dominan sob e los de conducción de calo y
los iscosos y se desa olla la con ección na u al, es e se denomina núme o de Rayleigh c í ico. El alo
eó ico de 𝑅𝑎𝑐=1708, que coincide muy bien con los alo es medidos, al ededo de 1708±50, pa a el caso
de dos placas in ini as.
Muchos in es igado es ealiza on mediciones de es e núme o de Rayleigh C í ico pa a ca idades ini as con
di e en es elaciones de aspec o y dis in as condiciones de con o no en las pa edes, po an o, exis e una g an
a iedad de bibliog a ía pa a compa a los esul ados.
Finalmen e comen a que ecien emen e se descub ió que los pa ones con ec i os obse ados po Béna d no
e an debidos a la di e encia de ensión supe icial que gene a la di e encia de empe a u as, no obs an e, es e
enómeno sigue denominándose comúnmen e como con ección de Rayleigh-Béna d.

17
2 ESTRUCTURA DEL TRABAJO
En el p esen e abajo se mues a un es udio del mé odo numé ico explíci o de MacCo mack
desa ollado po R.W.MacCo mack y su aplicación a di e sos p oblemas luidodinámicos pa a hace posible
su esolución median e un p og ama MATLAB ela i amen e sencillo.
En el capí ulo 3 analiza emos los undamen os de dicho mé odo, incluyendo es a egias que hacen que el
mé odo sea más es able. Se p esen a án las ecuaciones de Na ie -S okes en o ma conse a i a, es a
disposición de las ecuaciones es esencial pa a pode aplica es e esquema numé ico de o ma e icien e.
Además se incluye una pequeña aplicación de un caso unidimensional pa a ayuda a en ende como unciona
la mecánica del mé odo an es de comenza a esol e p oblemas más complejos y un ejemplo de un
p oblema bidimensional pa a comenza a en ende el mé odo aplicado a las ecuaciones de Na ie -S okes. Un
caso muy simple, una caja cuad ada llena de luido cuya apade a se desplaza con elocidad U cons an e. Se
simpli ican las ecuaciones de Na ie -S okes comple as mos adas al comienzo del apa ado y se aplican las
condiciones de con o no co espondien es. Es e ejemplo nos si e como e i icado del mé odo ya que ha
sido esuel o po muchos au o es y nos pe mi e compa a los esul ados.
El capí ulo 4, es una ex ensión del apa ado 4 del capí ulo 3, una ca idad ec angula ala gada, que se
puede asimila a un conduc o, con un cuad ado en medio que hace sus eces de obs áculo pa a el luido. Se
imponen las condiciones de con o no co espondien es al que el p oblema sea equi alen e al cuad ado
iajando a elocidad cons an e en el conduc o. El obje i o de es a aplicación se á el cálculo de la esis encia
gene ada po el cuad ado asi como la sus en ación si la hubiese. También amos a es udia la sensibilidad a
la elación de aspec o de la ca idad. Una ez lle ado a cabo es e análisis pod emos elegi una elación de
aspec o de la ca idad de al o ma que el esul ado se ap oxime a el cuad ado inme so en una co ien e
uni o me y en onces analiza emos el p oblema de ‘Vo ex shedding’ que comen amos an e io men e,
calculando la ecuencia con la que se emi en o bellinos aguas abajo.
En el capí ulo 5 se con empla la esolución de un p oblema pa ecido al an e io pe o en es e caso amos
a inclui el p oblema é mico. El cuad ado es a á a una empe a u a supe io a la del luido po an o el luido
ac ua á como e ige an e lle ándose el calo que pudiese desp ende el cuad ado. Se analiza án los pe iles
de empe a u a a dis in os núme os de P and l.
El capí ulo 6, ecoge un es udio del enómeno de ines abilidad é mica que se p oduce cuando una
ca idad es calen ada po su pa e in e io en p esencia de las ue zas másicas. Es a ines abilidad p oduce la
conocida como con ección na u al o con ección de Rayleigh-Béna d. Es e ejemplo es algo más complejo y
po an o se desa olla á un b e e a amien o eó ico del p oblema. Se ealiza án es udios en elación a las
dis in as celdas con ec i a que se o man al a ia la elación de aspec o y la in luencia de los núme os
adimensionales de nues o p oblema en la apa ición y igo osidad de la con ección ob enida. Es os análisis
se han ealizado sin hace uso de la ap oximación de Boussinesq, que es lo usual en la li e a u a. Es o
pe mi i á e el e ec o de la comp esibilidad en el égimen de con ección na u al.
18
3 MÉTODO DE MACCORMACK
3.1 Elección del mé odo
Exis en in inidad de mé odos numé icos pa a esol e p oblemas en mecánica de luidos, cada uno iene
sus en ajas y sus incon enien es. Según R.W. MacCo mack [28] odo p incipian e en CFD comienza
u ilizando mé odos de di e encias ini as po la amilia idad con las ecuaciones di e enciales y después
paula inamen e ienden a mé odos basados en olúmenes ini os. Como el obje i o p imo dial de es e
documen o es la aplicación de un mé odo numé ico a p oblemas de con ección pa a que pueda se
u ilizado po alumnos de asigna u as de mecánica de luidos hemos c eido mas ap opiado comenza con
un mé odo de di e encias ini as. Exis en muchos mé odos de di e encias ini as como los mé odos de
Lax-F ied ichs, Lax-Wend o o MacCo mack en e o os.
Hemos escogido el mé odo de MacCo mack po que además de p opo ciona esul ados más que
sa is ac o ios pa a di e sos p oblemas luidodinámicos, es el más sencillo de p og ama y po an o una
muy buena elección pa a inicia se en la mecánica de luidos compu acional.
También el Mé odo de MacCo mack admi e una ex ensión a p oblemas idimensionales bas an e sencilla
lo cual deja abie a una con inuación de es e abajo pa a esol e p oblemas luidodinámicos
idimensionales.
3.2 Fundamen oTeó ico
Es e mé odo consis e en dos e apas, una de p edicción y o a co ec o a. La p ime a consis e en una
ap oximación de las de i adas espaciales median e la di e encia p og esi a y la segunda median e la
di e encia eg esi a. Después se ealiza una media a i mé ica en e ambas gene ando una ap oximación
de segundo o den equi alen e a la di e encia cen al como se demos a á en el siguien e apa ado. Las
de i adas empo ales se ap oximan median e la di e encia p og esi a. Véase el siguien e esquema
gene al:
𝐸𝑐𝑢𝑎𝑐𝑖ó𝑛: 𝜕(𝑈)
𝜕𝑡 +𝜕(𝐹)
𝜕𝑥 +𝜕(𝐺)
𝜕𝑦 +𝑆=0
( 3-1)
𝑃𝑟𝑒𝑑𝑖𝑐𝑡𝑜𝑟:𝑈𝑖,𝑗
𝑛+1
=𝑈𝑖,𝑗
𝑛−𝛥𝑡(𝐹𝑖+1,𝑗
𝑛−𝐹𝑖,𝑗
𝑛
𝛥𝑥 )−𝛥𝑡(𝐺𝑖,𝑗+1
𝑛−𝐺𝑖,𝑗
𝑛
𝛥𝑦 )−𝛥𝑡𝑆
( 3-2)
𝐶𝑜𝑟𝑟𝑒𝑐𝑡𝑜𝑟:𝑈𝑖,𝑗
𝑛+1=12(𝑈𝑖,𝑗
𝑛+1
+𝑈𝑖,𝑗
𝑛−𝛥𝑡(𝐹𝑖,𝑗
𝑛+1
−𝐹𝑖−1,𝑗
𝑛+1

𝛥𝑥 )−𝛥𝑡(𝐺𝑖,𝑗
𝑛+1
−𝐺𝑖,𝑗−1
𝑛+1

𝛥𝑦 )−𝛥𝑡𝑆)
( 3-3)
Como se mues a en el esquema an e io la ecuación debe de ans o ma se en o ma conse a i a pa a la
ejecución de es e mé odo. Se puede in e cambia el o den de las ap oximaciones eg esi a y p og esi a
en el esquema supe io , es más, es ecomendable in e cambia los al e na i amen e pa a e i a di ecciones
espaciales p e e en es. En la siguien e abla se mues an las posibles combinaciones:
19
Tabla ( 3-1 )
P edic o
Co ec o
Di ección-x
Di ección-y
Di ección-x
Di ección-y
DP
DP
DR
DR
DB
DR
DP
DP
DP
DR
DR
DP
DR
DP
DP
DR
En muchos casos al e na es os ciclos de de i adas p opo ciona una mejo a en la es abilidad numé ica.
No exis e una demos ación pa a la con e gencia del mé odo en p oblemas bidimensionales, pe o la
siguien e condición se ha mos ado e ec i a en la p ác ica como e emos más adelan e:
𝛥𝑡≤1
|𝑢|
𝛥𝑥+|𝑣|
𝛥𝑦+𝑐√1
𝛥𝑥2+1
𝛥𝑦2
( 3-4)
Se conoce como condición de Cou an -F ied ich-Le y (CFL)
Es in e esan e no a que es e mé odo puede o igina un also es ado Pseudo- ansi o io al alcanza la unidad
de edondeo de la compu ado a y no e mina de con e ge debido a sus dos e apas.
Como se puede ap ecia en la ecuación an e io , pa a la aplicación de es e mé odo es esencial que las
ecuaciones es én en o ma conse a i a. A con inuación se mues an las ecuaciones de Na ie -S okes
comple as pa a el p oblema bidimensional en es a o ma:
𝜕(𝑈)
𝜕𝑡 +𝜕(𝐹)
𝜕𝑥 +𝜕(𝐺)
𝜕𝑦 +𝑆=0
( 3-5)
𝑈=[𝜌
𝜌𝑢
𝜌𝑣
𝐸]=[𝑈1
𝑈2
𝑈3
𝑈4]
( 3-6)
𝐹=
[
𝜌𝑢
𝜌𝑢2+𝑝−𝜏𝑥𝑥
𝜌𝑢𝑣−𝜏𝑥𝑦
(𝐸+𝑝)𝑢−𝑢𝜏𝑥𝑥−𝑣𝜏𝑥𝑦+𝑞𝑥
]
( 3-7)
20
𝐺=
[
𝜌𝑣
𝜌𝑢𝑣−𝜏𝑥𝑦
𝜌𝑣2+𝑝−𝜏𝑦𝑦
(𝐸+𝑝)𝑣−𝑢𝜏𝑥𝑦−𝑣𝜏𝑦𝑦+𝑞𝑦
]
( 3-8)
El ec o S, que no con iene de i adas espaciales ni empo ales, se denomina é mino uen e y se ese a
pa a las ue zas ex e nas que puedan apa ece en el luido como pueden se la g a edad o adición de calo .
Po ejemplo:
𝑆=[00
𝜌𝑔
𝜌𝑣𝑔]
( 3-9)
Una ez calculados los alo es de U se pueden ‘decodi ica ’ las a iables p ima ias de la siguien e o ma:
𝜌=𝑈1
𝑢=𝑈2
𝑈1
𝑣=𝑈3
𝑈1
( 3-10)
𝑒=𝑈4
𝑈1−(𝑈2
𝑈1)2+(𝑈3
𝑈1)2
2
( 3-11)
Que jun o a las ecuaciones de es ado nos pe mi i án de e mina odas las a iables luido-dinámicas. En es e
manual se conside a á el caso de un gas pe ec o:
𝑒=𝑐𝑣𝑇
𝑝=𝜌𝑅𝑔𝑇
( 3-12)
3.3 Ap oximación numé ica en mallas equiespacidas
El esquema de MacCo mack consis e en u iliza una e apa de di e encia p og esi a y o a e apa de
di e encia eg esi a y ealiza una combinación en e ambas. Es as es imaciones son cada una de p ime
o den pe o la combinación de ambas es equi alen e a una di e encia cen al que cómo e emos a
con inuación es de Segundo o den.
Las ap oximaciones más comunes pa a de i adas de p ime o den son:
21
 Di e encia p og esi a:
𝑑𝜙
𝑑𝑥≈𝜙𝑖+1−𝜙𝑖
Δ𝑥
( 3-13)
 Di e encia eg esi a:
𝑑𝜙
𝑑𝑥≈𝜙𝑖−𝜙𝑖−1
Δ𝑥
( 3-14)
 Di e encia cen al:
𝑑𝜙
𝑑𝑥≈𝜙𝑖+1−𝜙𝑖−1
2Δ𝑥
( 3-15)
Si se ealiza una expansión en se ie de Taylo de la unción 𝜙 al ededo del pun o 𝑥𝑖 ob end emos:
𝜙𝑖+1=𝜙𝑖+𝛥𝑥𝑑𝜙
𝑑𝑥𝑖+𝛥𝑥2
2𝑑2𝜙
𝑑𝑥2𝑖+𝛥𝑥3
6𝑑3𝜙
𝑑𝑥3𝑖+⋯
( 3-16)
𝜙𝑖=𝜙𝑖
( 3-17)
𝜙𝑖+1=𝜙𝑖−𝛥𝑥𝑑𝜙
𝑑𝑥𝑖+𝛥𝑥2
2𝑑2𝜙
𝑑𝑥2𝑖−𝛥𝑥3
6𝑑3𝜙
𝑑𝑥3𝑖+⋯
( 3-18)
Subs i uyendo es os desa ollos en las ap oximaciones an e io es enemos que:
 Di e encia p og esi a:
𝜙𝑖+1−𝜙𝑖
Δ𝑥 =𝑑𝜙
𝑑𝑥𝑖+Δ𝑥
2𝑑2𝜙
𝑑𝑥2𝑖+⋯
( 3-19)
El e o de uncamien oes de p ime o den.
 Di e encia eg esi a:
𝑑𝜙
𝑑𝑥=𝜙𝑖−𝜙𝑖−1
Δ𝑥
( 3-20)
𝜙𝑖−𝜙𝑖−1
Δ𝑥 =𝑑𝜙
𝑑𝑥𝑖−Δ𝑥
2𝑑2𝜙
𝑑𝑥2𝑖+⋯
( 3-21)
El e o de uncamien o es de p ime o den.

22
 Di e encia cen al:
𝜙𝑖+1−𝜙𝑖−1
2𝛥𝑥 =𝑑𝜙
𝑑𝑥𝑖−𝛥𝑥2
6𝑑3𝜙
𝑑𝑥3𝑖+⋯
( 3-22)
El e o de uncamien oes de Segundo o den.
Es in e esan e comen a que si sumamos la di e encia p og esi a y la di e encia eg esi a y di idios en e 2,
que es esencialmen e lo que hace el mé odo de MacCo mack, ob enemos la di e encia cen al que como se
acaba de comp oba es de segundo o den.
Pa a ac ualiza las condiciones de con o no en algunos casos necesi a emos aplica di e encias p og esi as y
eg esi as de segundo o den:
 Di e encia p og esi a:
𝑑𝜙
𝑑𝑥𝑖=−3𝜙𝑖−4𝜙𝑖+1+𝜙𝑖+2
2Δ𝑥 +𝑂(Δ𝑥2)
( 3-23)
 Di e encia eg esi a:
𝑑𝜙
𝑑𝑥𝑖=3𝜙𝑖−4𝜙𝑖−1+𝜙𝑖−2
2Δ𝑥 +𝑂(Δ𝑥2)
( 3-24)
En algunos casos es necesa io ealiza ap oximaciones de de i das de o den supe io , po ejemplo, pa a
ac ualiza condiciones de con o no o pa a compu a los é minos co espondien es a e ec os de iscosidad.
Las más ele an es son:
 Di e encia cen al:
𝜙𝑖+1−2𝜙𝑖+𝜙𝑖−1
Δ𝑥2
( 3-25)
 Di e encia p og esi a:
𝜙𝑖+2−2𝜙𝑖+1+𝜙𝑖
Δ𝑥2
( 3-26)
 Di e encia eg esi a:
𝜙𝑖−2𝜙𝑖−1+𝜙𝑖−2
Δ𝑥2
( 3-27)
Que con un p ocedmien ocompe amen e análogo al ealizado en el apa ado an e io se puede demos a que
la di e encia cen al es de segundo o den de p ecisión y las o as dos di e encias únicamen e de p ime o den.
En algunas ocasiones se á p eciso ap oxima de i adas de segundo o den con p ecisión de o den dos lo que
ealiza emos de la siguien e o ma:
23
 Di e encia p og esi a:
𝑑2𝜙
𝑑𝑥2𝑖=−2𝜙𝑖−5𝜙𝑖+1+4𝑢𝑖+2−𝜙𝑖+3
Δ𝑥2
( 3-28)
 Di e encia eg esi a:
𝑑2𝜙
𝑑𝑥2𝑖=2𝜙𝑖−5𝜙𝑖−1+4𝑢𝑖−2−𝜙𝑖−3
Δ𝑥2
( 3-29)
 Di e encia c uzada p og esi a:
𝜕2𝜙
𝜕𝑥𝜕𝑦𝑖,𝑗=−−(𝜙𝑖+2,𝑗−1−𝜙𝑖+2,𝑗+1)+4(𝜙𝑖+1,𝑗−1−𝜙𝑖+1,𝑗+1)−3(𝜙𝑖,𝑗−1−𝜙𝑖,𝑗+1)
4ΔxΔy +𝑂(Δ𝑥2,Δ𝑦2,Δ𝑥Δ𝑦)
( 3-30)
 Di e encia c uzada eg esi a:
𝜕2𝜙
𝜕𝑥𝜕𝑦𝑖,𝑗=−−(𝜙𝑖−2,𝑗+1−𝜙𝑖−2,𝑗−1)+4(𝜙𝑖−1,𝑗+1−𝜙𝑖−1,𝑗−1)−3(𝜙𝑖,𝑗+1−𝜙𝑖,𝑗−1)
4ΔxΔy +𝑂(Δ𝑥2,Δ𝑦2,Δ𝑥Δ𝑦)
( 3-31)
Con es as sencillas ó mulas pod emos aplica las condiciones de con o no y las ap oximaciones que
equie en las ecuaciones.
Las ecuaciones di e enciales que se mues an en el apa ado an e io han de se disc e izadas pa a pode
aplica es e mé odo de di e encias ini as. Se ha escogido una malla de (𝑁𝑥+1)×(𝑁𝑦+1) nodos donde
se e alua án las magni udes luidodinámicas necesa ias. La dis ancia en e cada nodo se denomina á Δ𝑥 en
di ección x y Δ𝑦 en di ección y.
Vamos a di e encia en e nodos in e io es y nodos on e a, los nodos in e io es (2:𝑁𝑥)×(2:𝑁𝑦) que se á
donde se aplica án y calcula an las magni udes luidodinámicas p o enien es de las ecuaciones de Na ie -
S okes y los nodos on e a donde se á necesa io impone condiciones de con o no, como e emos más
adelan e.
Es undamen al aplica las ecuaciones en o ma compac a ya que la condición CFL nos limi a mucho el
in e alo empo al máximo que podemos escoge gene ando unos iempos de compu ación
conside ablemen e ele ados. Po an o la inco po ación de bucles nos demo a ía aún mas lo cual p oduci ía
iempos de compu ación inabo dables.
A con inuación mos amos como se han de implemen a las ecuaciones, u iliza emos la ecuación de
con inuidad a modo de ejemplo.
𝜕𝜌
𝜕𝑡=−𝜕(𝜌𝑢)
𝜕𝑥 −𝜕𝜌𝑣
𝜕𝑦
( 3-32)
24
Un aspec o cla e que simpli ica mucho la implemen ación es la de inición de a iables auxilia es sob e las
cuales aplicamos di ec amen e la de i ada, po ejemplo, en la ecuación de con inuidad es as a iables se ían
𝜌𝑢 𝑦 𝜌𝑣 que pueden se calculadas i ialmen e a pa i de la elocidad y la densidad.
En p ime luga amos a conside a el é mino p edic o , en es a e apa u ilizábamos di e encia p og esi a.
Reco demos que an o la densidad como el es o de las a iables luidodinámicas dependen del iempo y de
la posición en el espacio, po an o los é minos de la ecuación queda ían como sigue.
𝜕𝜌
𝜕𝑡≅𝜌𝑖,𝑗
𝑡+1−𝜌𝑖,𝑗𝑡
Δ𝑡
( 3-33)
𝜕𝜌𝑢
𝜕𝑥≅𝜌𝑢𝑖+1,𝑗
𝑡−𝜌𝑢𝑖,𝑗𝑡
Δ𝑥
( 3-34)
𝜕𝜌𝑣
𝜕𝑦≅𝜌𝑣𝑖,𝑗+1
𝑡−𝜌𝑣𝑖,𝑗𝑡
Δ𝑦
( 3-35)
Quedando la ecuación:
𝜌𝑖,𝑗
𝑡+1−𝜌𝑖,𝑗𝑡
Δ𝑡 =−𝜌𝑢𝑖+1,𝑗
𝑡−𝜌𝑢𝑖,𝑗𝑡
Δ𝑥 −𝜌𝑣𝑖,𝑗+1
𝑡−𝜌𝑣𝑖,𝑗𝑡
Δ𝑦
( 3-36)
En la cual obse amos que odos los alo es son conocidos en el ins an e siendo la única incogni a la
densidad en el ins an e +1 que podemos calcula como:
𝜌𝑖,𝑗
𝑡+1=𝜌𝑖,𝑗𝑡−Δ𝑡𝜌𝑢𝑖+1,𝑗
𝑡−𝜌𝑢𝑖,𝑗𝑡
Δ𝑥 −Δ𝑡𝜌𝑣𝑖,𝑗+1
𝑡−𝜌𝑣𝑖,𝑗𝑡
Δ𝑦
( 3-37)
Que en o ma ma icial se puede exp esa como:
𝜌(2:𝑁𝑥,2:𝑁𝑦)𝑡+1
=𝜌(2:𝑁𝑥,2:𝑁𝑦)𝑡−Δ𝑡
Δ𝑥(𝜌𝑢(3:𝑁𝑥+1,2:𝑁𝑥)𝑡−𝜌𝑢(2:𝑁𝑥,2:𝑁𝑥)𝑡)
−Δ𝑡
Δ𝑦(𝜌𝑣(2:𝑁𝑥,3:𝑁𝑥+1)𝑡−𝜌𝑣(2:𝑁𝑥,2:𝑁𝑥)𝑡)
( 3-38)
A es e alo lo enomb a emos como 𝜌𝑖,𝑗
𝑡+1






po se el é mino p edic o . Aho a ealiza emos la misma
ope ación pe o pa a el é mino co ec o , en el que eco demos que se ap oximaban sus é minos median e la
di e encia eg esi a.
𝜕𝜌
𝜕𝑡≅𝜌𝑖,𝑗
𝑡+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 obse amos que odos los alo es son conocidos en el ins an e 𝑡+1







siendo la única incogni a la
densidad en el ins an e +1 que podemos calcula como:
𝜌𝑖,𝑗
𝑡+1=𝜌𝑖,𝑗𝑡+1






−Δ𝑡𝜌𝑢𝑖,𝑗
𝑡+1






−𝜌𝑢𝑖−1,𝑗𝑡+1






Δ𝑥 −Δ𝑡𝜌𝑣𝑖,𝑗
𝑡+1






−𝜌𝑢𝑖−1,𝑗𝑡+1






Δ𝑦
( 3-43)
A es e alo de la densidad lo enomb a emos como 𝜌𝑖,𝑗
𝑡+1∗ po se el co ec o .
𝜌𝑖,𝑗
𝑡+1=𝜌𝑖,𝑗
𝑡+1∗+𝜌𝑖,𝑗
𝑡
2
( 3-44)
Que en o ma ma icial se puede exp esa 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í únicamen e se ha mos ado el desa ollo de la ecuación de con inuidad po se común en odos los
ejemplos, en cada aplicación encon a án un código MATLAB con odas las ecuaciones implemen adas.
An es de comenza p oblemas más complejos es con enien e e un ejemplo unidimensional pa a cap a la
esencia del mé odo.
32
%====================ECUACIÓN DE ONDA====================%
%LAS VARIABLES P_ SE CORRESPONDEN A LOS TÉRMINOS DEL PREDICTOR%
%========================================================%
close all;clea 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
D =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===========================%
o =1:T/D
P_u(2:nx)=u(2:nx)-(D /Dx)*c*(u(3:nx+1)-u(2:nx));%%PREDICTOR
P_u(1)=u(1)-(D /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)+(D /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)-(D /Dx)*c*(P_u(2:nx)-P_u(1:nx-
1)))*0.5;%%CORRECTOR
u(1)=-sin( *D );%%C.C x=0
u(nx+1)=(u(nx+1)+P_u(nx+1)+(D /2/Dx)*c*(-P_u(nx-1)+4*P_u(nx)-
3*P_u(nx+1)))*0.5;%%C.C x=2PI
end
plo (linspace(0,2*pi,leng h(u)),u,'*')%%PLOT SOLUCIÓN NUMÉRICA
hold on
p=sin(- *D *c:2*pi/1000:2*pi- *D *c);%%SOLUCIÓN ANALÍTICA
plo ((linspace(0,2*pi,leng h(p))),p,' ');%%PLOT SOLUCIÓN ANALÍTICA
legend('Solución Numé ica','Solución Analí ica')
i le('N_x=50 T=5s')
xlabel('x')
ylabel('u')
%========================================================%

33
3.4 Aplicación: Ca idad cuad ada, con la apade a desplazándose a elocidad
cons an e.
Una ez analizado y comp endido el p oblema unidimensional amos a ex ende el mé odo a dos
dimensiones. El p oblema que se mues a a con inuación iene siendo usado como alidación de códigos de
CFD ya que posee una geome ía simple, las condiciones de con o no ambién lo son y lo más impo an e, ha
sido ep oducido po a ios au o es y exis e bibliog a ía con la que compa a .
Figu a( 3-7 )
A con inuación, se mues an las ecuaciones de Na ie -S okes en las que se han sup imido los é minos
co espondien es a las ue zas másicas debido a que es as son desp eciables en e a las ue zas de ine cia
pa a los casos analizados:
𝐹𝑟2=𝑈2
𝑔𝐷≫1
( 3-65)
También amos a conside a que como el luido se encuen a inicialmen e a empe a u a cons an e y odas
las pa edes se encuen an a esa empe a u a, las a iaciones de empe a u a se án muy pequeñas si las
elocidades son pequeñas, como e emos a con inuación pudiendo conside a el caso iso e mo, ealizando
es imaciones de ó denes de magni ud en la ecuación de la en alpía:
𝜌𝑐𝑝𝜕𝑇
𝜕𝑡+𝜌𝑐𝑝𝑣
󰇍
󰇍
.∇𝑇=𝜕𝑝
𝜕𝑡+𝑣
󰇍
󰇍
.∇𝑝+∇.(𝜅∇T)+τ′:∇𝑣
( 3-66)
Si se in oducen los siguien es núme os adimensionales:
𝑅𝑒=𝜌𝑈𝐷
𝜇,𝑃𝑟=𝜈𝛼, 𝑀= 𝑢
𝛾𝑅𝑔𝑇0,
( 3-67)
34
Se ob iene:
 Inc emen os de empe a u a debido a a iaciones de p esión.
Δ𝑇𝑝
𝑇0~𝑈2
𝐶𝑝𝑇0~𝑀2
( 3-68)
Donde se ha enido en cuen a que: Δ𝑝 ~𝑈2
𝐷
 Inc emen os de empe a u a debido a conducción de calo .
Δ𝑇𝑞
𝑇0~𝑇1−𝑇2
𝑇0𝑃𝑟𝑅𝑒
( 3-69)
 Inc emen os de empe a u a debido a disipación iscosa.
Δ𝑇𝑑𝑖𝑠
𝑇0~𝑈2
𝐶𝑝𝑇0𝑅𝑒~𝑀2
𝑅𝑒
( 3-70)
En es as es exp esiones ap eciamos como si el núme o de Mach es pequeño, las empe a u as de las pa edes
son ap oximadamen e iguales y el núme o de Reynolds no es pequeño, la ap oximación po el caso iso e mo
es muy adecuada.
Po an o, las ecuaciones de Na ie -S okes bidimensionales se pueden exp esa como:
𝐶𝑜𝑛𝑡𝑖𝑛𝑢𝑖𝑑𝑎𝑑: 𝜕𝜌
𝜕𝑡+𝜕(𝜌𝑢)
𝜕𝑥 +𝜕(𝜌𝑣)
𝜕𝑦 =0
( 3-71)
𝐶𝑎𝑛𝑡𝑖𝑑𝑎𝑑 𝑑𝑒 𝑚𝑜𝑣𝑖𝑚𝑖𝑒𝑛𝑡𝑜 𝑥:𝜕(𝜌𝑢)
𝜕𝑡 +𝜕(𝜌𝑢2+𝑝−4𝜇
3𝑑𝑢
𝑑𝑥−𝜇3𝑑𝑣
𝑑𝑦)
𝜕𝑥 +𝜕(𝜌𝑢𝑣−𝜇𝑑𝑢
𝑑𝑦)
𝜕𝑦 =0
( 3-72)
𝐶𝑎𝑛𝑡𝑖𝑑𝑎𝑑 𝑑𝑒 𝑚𝑜𝑣𝑖𝑚𝑖𝑒𝑛𝑡𝑜 𝑦:𝜕(𝜌𝑣)
𝜕𝑡 +𝜕(𝜌𝑢𝑣−𝜇𝑑𝑣
𝑑𝑥)
𝜕𝑥 +𝜕(𝜌𝑣2+𝑝−4𝜇
3𝑑𝑣
𝑑𝑦−𝜇3𝑑𝑢
𝑑𝑥)
𝜕𝑦 =0
( 3-73)
ecuación de es ado: p=𝑅𝑔𝑇0𝜌
( 3-74)
Siendo c la elocidad del sonido en el medio.
35
Como hemos is o an e io men e a bajas elocidades la empe a u a debe ía se ap oximadamen e cons an e
po an o es as ecuaciones ap oximan el lími e iso e mo.
Es con enien e in oduci las siguien es a iables adimensionales:
𝑥∗=𝑥𝐷,𝑦∗=𝑦𝐷,𝑢∗=𝑢𝑈,𝑣∗=𝑣𝑈,𝑡∗=𝑈
𝐷𝑡,𝜌∗=𝜌
𝜌0,𝑝∗=𝑝
𝜌0𝑈2
( 3-75)
Siendo D un lado de la ca idad, U la elocidad de la apa y 𝜌0 una densidad de e e encia.
A pa i de aho a se sup imi á el as e isco de las a iables adimensionales po simplicidad. Las ecuaciones
esul an es 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)
Suje o a las condiciones de con o no:
𝑢(0,𝑦)=0 𝑣(0,𝑦)=0 𝑢(1,𝑦)=0 𝑣(1,𝑦)=0
( 3-80)
𝑢(𝑥,0)=0 𝑣(𝑥,0)=0 𝑢(𝑥,1)=1 𝑣(𝑥,1)=0
( 3-81)
Un ac o impo an e a la ho a de implemen a sa is ac o iamen e el mé odo es impone co ec amen e las
condiciones de con o no en densidad o equi alen emen e en p esión. Vamos a ealiza es o u ilizando la
ecuación de con inuidad:
𝑇𝑎𝑝𝑎 𝑖𝑛𝑓𝑒𝑟𝑖𝑜𝑟: 𝜕𝜌
𝜕𝑡=−𝜕(𝜌𝑣)
𝜕𝑦
( 3-82)
36
𝑇𝑎𝑝𝑎 𝑠𝑢𝑝𝑒𝑟𝑖𝑜𝑟: 𝜕𝜌
𝜕𝑡=−𝜕𝜌
𝜕𝑥−𝜕(𝜌𝑣)
𝜕𝑦
( 3-83)
𝑇𝑎𝑝𝑎 𝑖𝑧𝑞𝑢𝑖𝑒𝑟𝑑𝑎: 𝜕𝜌
𝜕𝑡=−𝜕(𝜌𝑢)
𝜕𝑥
( 3-84)
𝑇𝑎𝑝𝑎 𝑑𝑒𝑟𝑒𝑐ℎ𝑎: 𝜕𝜌
𝜕𝑡=−𝜕(𝜌𝑢)
𝜕𝑥
( 3-85)
Las ecuaciones de Na ie S okes se encuen an disc e izadas en el apa ado 2, únicamen e end emos que
desp ecia las ue zas másicas y al conside a el caso iso e mo pod emos ob ia la ecuación de la ene gía.
Vamos a analiza en más de alle la disc e ización de las condiciones de con o no en densidad que son algo
más complejas.
A pa i de es as exp esiones podemos calcula la densidad en el con o no de la caja median e las siguien es
ap oximaciones numé icas que imos en el segundo apa ado. Pa a la de i ada empo al pod emos u iliza la
di e encia p og esi a.
𝜕𝜌
𝜕𝑡≅𝜌𝑖,𝑗
𝑡+1−𝜌𝑖,𝑗𝑡
Δ𝑡
( 3-86)
Sin emba go no se puede u iliza la di e encia p og esi a/ eg esi a indis in amen e pa a una misma e apa ya
que al es a en la on e a es as di e encias eque i ían nodos ue a del dominio compu acional conside ado,
en o as palab as excede íamos la dimensión de la ma iz. Pues o que no se puede segui el guión gene al del
mé odo, el cual le con ie e una p ecisión de segundo o den en el iempo y en el espacio, y que emos
man ene dicha p ecisión debemos aplica ap oximaciones de segundo o den. La di e encia cen al no se
puede u iliza ya que excede íamos nue amen e la dimensión de la ma iz po an o elegi emos la di e encia
p og esi a o eg esi a de al o ma que siemp e nos man engamos den o del dominio compu acional. Los
alo es de la densidad en las esquinas se co esponde án con los alo es ex emos de las apas la e ales lo
cual nos pe mi e ap oxima la de i ada espacial de la densidad que apa ece en la apa supe io median e la
di e encia cen al.
𝑇𝑎𝑝𝑎 𝑖𝑧𝑞𝑢𝑖𝑒𝑟𝑑𝑎: 𝜕𝜌𝑢
𝜕𝑥𝑖,𝑗≅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 es as exp esiones no es di ícil calcula los alo es de la densidad en el ins an e
𝑡+1







(𝑐𝑜𝑟𝑟𝑒𝑠𝑝𝑜𝑛𝑑𝑖𝑒𝑛𝑡𝑒 𝑎𝑙 𝑝𝑟𝑒𝑑𝑖𝑐𝑡𝑜𝑟) en la on e a 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 ins an e +1(co ec o ) 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 MacCo mack:
1. E apa de p edicción (Di e encia P og esi a)
En es a ase calculamos las a iables 𝜌𝑡+1






,𝜌𝑢𝑡+1






𝑦 𝜌𝑣𝑡+1






median e las ecuaciones
disc e izadas de con inuidad y de can idad de mo imien o.
2. Decodi icación de a iables
Se calcula u y a pa i de las a iables que ob enemos de las ecuaciones.
3. Condiciones de con o no.
Aho a aplica emos las condiciones de con o no en elocidad y en densidad, como se acaba de
explica .
4. Ac ualización de a iables
Es impo an e ac ualiza nues as a iables auxilia es
𝜌𝑢𝑡+1






,𝜌𝑣𝑡+1






,𝜌𝑢2𝑡+1






,𝜌𝑣2𝑡+1






𝑦 𝜌𝑢𝑣𝑡+1






despues de aplica las condiciones de con o no.
5. E apa de co ección (Di e encia eg esi a).
En es a ase calculamos las a iables 𝜌𝑡+1,𝜌𝑢𝑡+1 𝑦 𝜌𝑣𝑡+1 median e las ecuaciones
disc e izadas de con inuidad y de can idad de mo imien o.
6. Decodi icación de a iables
Se calcula u y a pa i de las a iables que ob enemos de las ecuaciones.
7. Condiciones de con o no.
Las condiciones de con o no de es a e apa se ealizan siguiendo el esquema an e io
8. Ac ualización de a iables
Es impo an e ac ualiza nues as a iables auxilia es
𝜌𝑢𝑡+1 ,𝜌𝑣𝑡+1,𝜌𝑢2𝑡+1 ,𝜌𝑣2𝑡+1 𝑦 𝜌𝑢𝑣𝑡+1 despues de aplica las condiciones de con o no.
9. Repe imos el p oceso has a alcanza el iempo deseado.
Un aspec o que aun no hemos analizado es el inc emen o de iempo que amos a selecciona . Hemos
comen ado en el p ime apa ado que, a pesa de no exis i una demos ación de una condición de es abilidad
pa a p oblemas mul idimensionales, la condición de Cou an -F ied ich-Le y unciona conside ablemen e
bien en la p ác ica, la eco damos:
𝛥𝑡≤1
|𝑢|
𝛥𝑥+|𝑣|
𝛥𝑦+𝑐√1
𝛥𝑥2+1
𝛥𝑦2
(
3-100)
39
Incluso mejo , pa a p oblemas que concie nen a la esolución de las ecuaciones de Na ie -S okes, Tannehill
e al., nos p opone la siguien e ó mula empí ica.
𝛥𝑡∗≤𝜎Δ𝑡
1+2/𝑅𝑒Δ
( 3-101)
Donde 𝑅𝑒Δ es el núme o de Reynolds basado en el mínimo de Δ𝑥 𝑦 Δ𝑦 y 𝜎=0.9, es un ac o de
segu idad. Obse amos que pa a alo es ele ados del Núme o de Reynolds basado en Δ las condiciones son
simila es, siendo la segunda algo más es ic i a.
Vol iendo a la condición CFL obse amos que la podemos simpli ica si enemos en cuen a que Δ𝑥 ha sido
elegido igual a Δ𝑦 y que el alo máximo de u y se á del mismo o den de magni ud.
𝛥𝑡≤Δ𝑥
2|𝑢|+𝑐√2
( 3-102)
Que obse amos que es más exigen e que la que ob eníamos en una a iable.
𝛥𝑡≤Δ𝑥
|𝑢|+𝑐
( 3-103)
A modo de ano ación y po si se quisiese op imiza el código R.W. MacCo mack hizo una e isión de su
mé odo en el que di ide su esquema en e apas unidimensionales, lo que como acabamos de e educe las
es icciones de iempo y mejo a aún más si com emplamos mallas en las que Δ𝑥 y Δ𝑦 sean muy di e en es
ya que se pod ía a anza en cada di ección con el inc emen o de iempo máximo pe mi ido. No obs an e,
amos a aplica el esquema o iginal po su simplicidad y buenos esul ados.
Vamos a analiza los casos de 𝑅𝑒=100 𝑦 𝑅𝑒=400 pa a 𝑀=0.05 pa a man ene nos siemp e en el
lími e incomp esible.
40
Figu a( 3-8 )
Figu a( 3-9 )
Podemos obse a como la apa ‘a as a’ al luido gene ando la o ación y o mando el emolino an
ca ac e ís ico. Las líneas de co ien e más oscu as co esponden a los esul ados de Hou y el campo ec o ial
de elocidades g is cla o al ob enido po Kundu, Cohen & Dowling.
41
Figu a( 3-10 )
En el de alle de la imagen se ap ecia como el cen o de la celda con ec i a se encuen a en [0,62 0,74] que
son muy p óximos a los alo es ob enidos po Hou e al.(1995) , ealizados median e el mé odo
la iceBol zman, que aco daban que el cen o se encon aba en [0,6196 0,7373].
Figu a( 3-11 )
48
𝑒𝑐𝑢𝑎𝑐𝑖ó𝑛 𝑑𝑒 𝑒𝑠𝑡𝑎𝑑𝑜: 𝑝= 𝜌
𝛾𝑀2
( 4-4)
En las p óximas aplicaciones po simplicidad se esc ibi án las ecuaciones de Na ie -S okes en no ación
ma icial.
El p oblema es á suje o a las condiciones de con o no:
𝑢(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 coo denada x del inicio del cuad ado, 𝑥𝑐𝑓 la coo denada x del inal del cuad ado y pa a la
coo denada y la nomencla u a es análoga.
Es necesa io impone condiciones de con o no en densidad (o p esión). Vamos a ealiza es o u ilizando la
ecuación de con inuidad:
𝑇𝑎𝑝𝑎 𝑖𝑛𝑓𝑒𝑟𝑖𝑜𝑟(𝑑(𝜌𝑢)
𝑑𝑥 =𝑑𝜌
𝑑𝑥)𝑑𝜌
𝑑𝑡=−𝑑𝜌
𝑑𝑥−𝑑(𝜌𝑣)
𝑑𝑦
( 4-9)
𝑇𝑎𝑝𝑎 𝑠𝑢𝑝𝑒𝑟𝑖𝑜𝑟(𝑑(𝜌𝑢)
𝑑𝑥 =𝑑𝜌
𝑑𝑥)𝑑𝜌
𝑑𝑡=−𝑑𝜌
𝑑𝑥−𝑑(𝜌𝑣)
𝑑𝑦
( 4-10)
𝑇𝑎𝑝𝑎 𝑖𝑧𝑞𝑢𝑖𝑒𝑟𝑑𝑎(𝑑(𝜌𝑣)
𝑑𝑦 =0)𝑑𝜌
𝑑𝑡=−𝑑(𝜌𝑢)
𝑑𝑥
( 4-11)
𝑇𝑎𝑝𝑎 𝑑𝑒𝑟𝑒𝑐ℎ𝑎(𝑑(𝜌𝑣)
𝑑𝑦 =0)𝑑𝜌
𝑑𝑡=−𝑑(𝜌𝑢)
𝑑𝑥
( 4-12)
Que obse amos que si aplicamos las di e encias p og esi as, eg esi as y cen ales ob enidas en el apa ado
3.2 pod emos calcula las densidades en la on e a de nues o canal.
Pa a las condiciones de densidad en el cuad ado, se puede ac ualiza o bien con la ecuación de con inuidad:
𝐿𝑎𝑑𝑜 𝑖𝑛𝑓𝑒𝑟𝑖𝑜𝑟(𝑑(𝜌𝑢)
𝑑𝑥 =0)𝑑𝜌
𝑑𝑡=−𝑑(𝜌𝑣)
𝑑𝑦
( 4-13)
𝐿𝑎𝑑𝑜 𝑠𝑢𝑝𝑒𝑟𝑖𝑜𝑟(𝑑(𝜌𝑢)
𝑑𝑥 =0)𝑑𝜌
𝑑𝑡=−𝑑(𝜌𝑣)
𝑑𝑦
( 4-14)

49
𝐿𝑎𝑑𝑜 𝑖𝑧𝑞𝑢𝑖𝑒𝑟𝑑𝑜(𝑑(𝜌𝑣)
𝑑𝑦 =0)𝑑𝜌
𝑑𝑡=−𝑑(𝜌𝑢)
𝑑𝑥
( 4-15)
𝑇𝑎𝑝𝑎 𝑑𝑒𝑟𝑒𝑐ℎ𝑜(𝑑(𝜌𝑣)
𝑑𝑦 =0)𝑑𝜌
𝑑𝑡=−𝑑(𝜌𝑢)
𝑑𝑥
( 4-16)
O bien desde la ecuación de can idad de mo imien o, que es más compleja, pe o en es e caso pa ece da
mejo es esul ados:
𝐿𝑎𝑑𝑜 𝑖𝑛𝑓𝑒𝑟𝑖𝑜𝑟 𝑑𝜌
𝑑𝑦=𝑀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 ap oximaciones es ablecidas en el capí ulo 3 se pueden despeja las a iables en el con o no
del cuad ado.
Pa a calcula los coe icien es de esis encia y sus en ación amos a u iliza las siguien es exp esiones:
𝐶𝑙=𝑙
12𝜌𝑈2𝐷=∫(𝑃𝑖𝑛𝑓−𝑃𝑠𝑢𝑝)𝑑𝑥+∫(𝑓𝑦𝑖𝑧𝑞+𝑓𝑦𝑑𝑒𝑟)𝑑𝑦+∫(𝑓𝑦𝑠𝑢𝑝+𝑓𝑦𝑖𝑛𝑓)𝑑𝑥
𝑖2
𝑖1
𝑗2
𝑗1
𝑖2
𝑖112𝜌𝑈2𝐷
( 4-21)
50
Siendo 𝑓𝑦 los es ue zos angenciales que ob enemos en las pa edes:
𝑓𝑦𝑖𝑧𝑞=−𝜇𝑑𝑣
𝑑𝑥𝑥=𝑖1
( 4-22)
𝑓𝑦𝑑𝑒𝑟=𝜇𝑑𝑣
𝑑𝑥𝑥=𝑖2
( 4-23)
𝑓𝑦𝑑𝑒𝑟=𝜇𝑑𝑣
𝑑𝑥𝑥=𝑖2
( 4-24)
𝑓𝑦𝑠𝑢𝑝=2𝜇𝑑𝑣
𝑑𝑦𝑦=𝑗2
( 4-25)
𝑓𝑦𝑖𝑛𝑓=−2𝜇𝑑𝑣
𝑑𝑦𝑦=𝑗1
( 4-26)
Con i iendo a nues as a iables adimensionales:
𝐶𝑙
=𝜌𝑈2𝐷∫(𝑃𝑖𝑛𝑓
∗−𝑃𝑠𝑢𝑝
∗)𝑑𝑥∗+𝜇𝑈∫(𝑓𝑦∗𝑖𝑧𝑞+𝑓𝑦∗𝑑𝑒𝑟)𝑑𝑦∗
𝑗2
𝐷
𝑗1
𝐷+𝜇𝑈∫(𝑓𝑦∗𝑠𝑢𝑝+𝑓𝑦∗𝑖𝑛𝑓)𝑑𝑥∗
𝑗2
𝐷
𝑗1
𝐷
𝑖2
𝐷
𝑖1
𝐷12𝜌𝑈2𝐷
=2∫(𝑃𝑖𝑛𝑓
∗−𝑃𝑠𝑢𝑝
∗)𝑑𝑥∗+2
𝑅𝑒∫ (𝑓𝑦∗𝑖𝑧𝑞+𝑓𝑦∗𝑑𝑒𝑟)𝑑𝑦∗+2
𝑅𝑒∫ (𝑓𝑦∗𝑠𝑢𝑝+𝑓𝑦∗𝑖𝑛𝑓)𝑑𝑥∗
𝑖2
𝐷
𝑖1
𝐷
𝑗2
𝐷
𝑗1
𝐷
𝑖2
𝐷
𝑖1
𝐷
( 4-27)
Lo cual iene sen ido ísico ya que un aumen o de la di e encia de p esiones en e las ca as supe io e in e io
con ibuye posi i amen e a la sus en ación. También es cohe en e que aumen a el núme o de Reynolds
p oduzca una pe dida de impo ancia de las ue zas de iscosidad.
Si enemos en cuen a que 𝑃∗=𝜌∗
𝑀2:
𝐶𝑙=2
𝑀2∫(𝜌𝑖𝑛𝑓
∗−𝜌𝑠𝑢𝑝
∗)𝑑𝑥∗+
𝑖2
𝐷
𝑖1
𝐷2
𝑅𝑒∫ (𝑓𝑦∗𝑖𝑧𝑞+𝑓𝑦∗𝑑𝑒𝑟)𝑑𝑦∗+2
𝑅𝑒∫ (𝑓𝑦∗𝑠𝑢𝑝+𝑓𝑦∗𝑖𝑛𝑓)𝑑𝑥∗
𝑖2
𝐷
𝑖1
𝐷
𝑗2
𝐷
𝑗1
𝐷
( 4-28)
Análogamen e pa a el coe icien e de esis encia:
51
𝐶𝑑=𝑑
12𝜌𝑈2𝐷=∫(𝑃𝑖𝑧𝑞−𝑃𝑑𝑒𝑟)𝑑𝑦+∫(𝑓𝑥𝑖𝑧𝑞+𝑓𝑥𝑑𝑒𝑟)𝑑𝑦+∫(𝑓𝑥𝑠𝑢𝑝+𝑓𝑥𝑖𝑛𝑓)𝑑𝑥
𝑖2
𝑖1
𝑗2
𝑗1
𝑗2
𝑗112𝜌𝑈2𝐷
( 4-29)
Siendo 𝑓𝑥 los es ue zos angenciales que ob enemos en las pa edes:
𝑓𝑥𝑖𝑧𝑞=−2𝜇𝑑𝑢
𝑑𝑥𝑥=𝑖1
( 4-30)
𝑓𝑥𝑑𝑒𝑟=2𝜇𝑑𝑢
𝑑𝑥𝑥=𝑖2
( 4-31)
𝑓𝑥𝑠𝑢𝑝=𝜇𝑑𝑢
𝑑𝑦𝑦=𝑗2
( 4-32)
𝑓𝑥𝑖𝑛𝑓=−𝜇𝑑𝑢
𝑑𝑦𝑦=𝑗1
( 4-33)
Con i iendo a nues as a iables adimensionales:
𝐶𝑑
=𝜌𝑈2𝐷∫(𝑃𝑖𝑧𝑞
∗−𝑃𝑑𝑒𝑟
∗)𝑑𝑦∗+𝜇𝑈∫(𝑓𝑥∗𝑖𝑧𝑞+𝑓𝑥∗𝑑𝑒𝑟)𝑑𝑦∗
𝑗2
𝐷
𝑗1
𝐷+𝜇𝑈∫(𝑓𝑥∗𝑠𝑢𝑝+𝑓𝑥∗𝑖𝑛𝑓)𝑑𝑥∗
𝑗2
𝐷
𝑗1
𝐷
𝑗2
𝐷
𝑗1
𝐷12𝜌𝑈2𝐷
=2∫ (𝑃𝑖𝑧𝑞
∗−𝑃𝑑𝑒𝑟
∗)𝑑𝑦∗+2
𝑅𝑒∫ (𝑓𝑥∗𝑖𝑧𝑞+𝑓𝑥∗𝑑𝑒𝑟)𝑑𝑦∗+2
𝑅𝑒∫ (𝑓𝑥∗𝑠𝑢𝑝+𝑓𝑥∗𝑖𝑛𝑓)𝑑𝑥∗
𝑖2
𝐷
𝑖1
𝐷
𝑗2
𝐷
𝑗1
𝐷
𝑗2
𝐷
𝑗1
𝐷
( 4-34)
Lo cual iene sen ido ísico ya que un aumen o de la di e encia de p esiones en e las ca as izquie da y
de echa con ibuye posi i amen e a la esis encia. También es cohe en e que aumen a el núme o de
Reynolds p oduzca una pe dida de impo ancia de las ue zas de iscosidad.
Si enemos en cuen a que 𝑃∗=𝜌∗
𝑀2:
𝐶𝑑=2
𝑀2∫(𝜌𝑖𝑧𝑞
∗−𝜌𝑑𝑒𝑟
∗)𝑑𝑦∗+
𝑗2
𝐷
𝑗1
𝐷2
𝑅𝑒∫ (𝑓𝑥∗𝑖𝑧𝑞+𝑓𝑥∗𝑑𝑒𝑟)𝑑𝑦∗+2
𝑅𝑒∫ (𝑓𝑥∗𝑠𝑢𝑝+𝑓𝑥∗𝑖𝑛𝑓)𝑑𝑥∗
𝑖2
𝐷
𝑖1
𝐷
𝑗2
𝐷
𝑗1
𝐷
( 4-35)
Comenza emos analizando el lujo al ededo del cuad ado pa a 𝑅𝑒=20 𝑦 𝑀=0.05
An es de analiza los coe icien es de esis encia y sus en ación amos a analiza el pa ón de lujo que se
ob iene en o no a la caja y que le ocu e a la co ien e al ebo dea a es a.
52
Figu a( 4-2 )
Figu a( 4-3 )
Obse amos como inmedia amen e después de bo dea la caja se o man dos o bellinos que pe manecen
adhe idos, y obse amos como en ningún caso se p oduce oscilación en la es ela, ob eniendo un égimen
es aciona io.
Vamos a ep esen a aho a los coe icien es de esis encia y sus en ación.
53
Figu a( 4-4 )
Obse amos como el alo no es idén icamen e nulo, es o se debe a asíme ias en la malla, no obs an e, se
puede obse a que si e inamos la malla es e alo a disminuyendo has a ap oxima se mucho a 0.
Figu a( 4-5 )
Pe cibimos con es a imagen como es as asime ías se han educido conside ablemen e has a un alo de 0.03,
ya muy p óximo a sus en ación nula que e a el compo amien o espe ado.

54
Figu a( 4-6 )
En la igu a 4-6 obse amos como el alo del Cd no es esidual como ocu ia con el Cl en es e caso como
e a de espe a exis e una esis encia impo an e al mo imien o ya que la co ien e se desp ende, p o ocando
las bu bujas de eci culación en so a en o p oduciendo una impo an e esis encia.
O o aspec o signi ica i o es el hecho de que pa a es e núme o de Reynolds apa ezca un égimen
es aciona io, algo que no se á e dad pa a cualquie núme o de Reynolds como e emos a con inuación.
Re=100.
Figu a( 4-7 )
55
Figu a( 4-8 )
Figu a( 4-9 )
En la secuencia de imágenes an e io ,o denadas empo almen e, se mues a como el luido comienza a o a
al ebo dea la caja o mándose un o bellino que se a con ec ando aguas abajo haciendo que la es ela
oscile lo que como e emos a con inuación p o oca á a iaciones empo ales en los coe icien es de
esis encia y de sus en ación.
56
Figu a( 4-10 )
Como comen ábamos an es, ob enemos una oscilación en el coe icien e de sus en ación, es o se debe a la
oscilación de la es ela como se ha mos ado an e io men e, aunque si que obse amos que debido a la
sime ía del p oblema es a oscilación es simé ica espec o a la ho izon al de una ampli ud ap oximadamen e
0.8, muy simila a la ob enida po Kundu, Cohen &Dowling, que ob enían un alo de 0.77.
Figu a( 4-11 )
57
Obse amos que como e a de espe a la oscilación de la es ela p oduce una a iación de la esis encia, no
obs an e, es a es mucho más educida en ampli ud y de ecuencia mucho más al a, obse amos como su
alo medio es e en o no a los 3.4, en la siguien e imagen se ap ecia con más cla idad. De nue o un alo
muy simila al p edicho po Kundu, Cohen &Dowling que ob enían un alo medio de 3.35.
Figu a( 4-12 )
Es in e esan e compa a es e alo medio ob enido con el del apa ado an e io y obse a que se ha educido
p ác icamen e a la mi ad, lo cual iene sen ido ísico ya que al aumen a el núme o de Reynolds disminuye la
impo ancia ela i a de las ue zas de iscosidad y po an o de la esis encia.
Vamos a ealiza un análisis de sensibilidad de la elación de aspec o de la ca idad espec o del cuad ado
pa a un alo del núme o de Reynolds de 20 y después de un iempo adimensional de 15, cuando ya se ha
alcanzado el égimen es aciona io, en elación al coe icien e de esis encia. La posición ela i a del cuad ado
se ha man enido cons an e.
.
Tabla ( 4-1 )
𝑳/𝑫
𝑯/𝑫
𝑪𝒅
35
3
5.5703
100
3
6.0650
La condición de con o no que se ha pues o en la salida del conduc o es de elocidad uni o me igual a
elocidad de en ada, po an o, al aumen a la longi ud de la ca idad se iende al caso eal ya que los e ec os
pe u bado es p oduc o de ebo dea el obje o ya son desp eciables. En cambio si educimos la longi ud de la
ca idad obligamos a los e ec os pe u bado es de la es ela a amo igua se an es de lo que lo ha ían en un caso
eal.
64
Figu a( 4-22 )
Obse amos como en la igu a 4-15 co espondien e a 𝑅𝑒=2 la co ien e pe manece adhe ida, mien as
que en la igu a 4-16 co espondien e a 𝑅𝑒=5 es o ya no ocu e y se comienzan a o ma o bellinos que se
man ienen adhe idos a so a en o, sin p o oca la oscilación de la es ela. En la igu a 4-17 ya se obse an los
o bellinos de mayo in ensidad pe o se man ienen adhe idos. En la igu a 4-18 ya obse amos como se
comienza a ompe la sime ía y se comienza a emi i o icidad aguas abajo, el alo del núme o de
Reynolds c í ico pa a el cual comienza la emisión de o bellinos se á muy p óximo pe o lige amen e in e i o
a 𝑅𝑒=40 . Finalmen e en la igu a 4-19 se ap ecia muy cla amen e como se es án emi iendo o bellinos
aguas abajo p o ocando la oscilación de la es ela que an es comen ábamos.
Obse amos como los esul ados son muy pa ejos a los del cilind o ci cula in ini o, p obablemen e en
nues o caso los umb ales sean lige amen e in e io es, debido a que los angulos ec os en los que se unen las
ca as de nues o cilind o acili an el desp endimien o de la co ien e.
Llegados a es e pun o nos su ge la siguien e p eocupación, ¿Con que ecuencia se desp enden o bellinos
pa a dis in os núme os de Reynolds? Como hemos is o an e io men e es o se puede es udia median e el
núme o de S ouhal que es equi alen e a una ecuencia adimensional. La adimensionalización que hemos
ealizado no puede se más adecuada pa a el es udio de es e enómeno ya que pa a calcula la ecuencia
adimensional y po ende el núme o de S ouhal únicamen e debemos calcula el pe iodo adimensional he
in e i lo, es o se puede hace po ejemplo obse ando las oscilaciones del coe icien e de sus en ación y
ob eniendo su pe iodo.
En la igu a 4-22 obse amos los esul ados ob enidos median e el esquema de MacCo mack pa a una
ca idad de elaciones de aspec o 𝐿𝐷=35 y 𝐻
𝐷=6 los cuales debe ían ap oxima conside ablemen e bien a
un cilind o inme so en una co ien e uni o me jun o con los esul ados de la bibliog a ía.

65
Figu a( 4-23 )
Donde se obse a como el acue do con los da os ob enidos de di e osos au o es es muy ele ado. Pa a
ob ene una p ecisión incluso mayo podíamos amplia el dominio compu acional e impone las condiciones
de con o no más lejos del cilind o. Hemos op ado po es a solución de comp omiso po que p opo cionaba
esul ados que se ajus aban a la bibliog a ía sin eque i un iempo de compu ación excesi amen e ele ado.
Como hemos comen ado en la in oducción en la ingenie ía la a iación de sus en ación puede p oduci
p oblemas es uc u ales en an enas de ehículos, o es de e ige ación, ascacielos, en de ini i a es uc u as
que ean una co ien e ap oximadamen e uni o me. Pa a e i a la o u a po a iga de es os elemen os suelen
inclui ale as pa a a a de educi es as oscilaciones ya que si se e i e que los o bellinos es én en con ac o
di ec o la ines abilidad se educe conside ablemen e.
A con inuación mos amos dos imágenes una co espondien e al lujo al ededo de un cilind o ci cula , y en
la siguien e el cilind o inco po a una ale a pa a e i a la o mación de ó ices.
Figu a( 4-24 )
Se ap ecia cla amen e la educción en la oscilación de la es ela. Vamos a mos a el e ec o pa a nues o
cuad ado.
66
Figu a( 4-25 )
Figu a( 4-26 )
67
Figu a( 4-27 )
Se ap ecia cla amen e como la inclusión de la ale a p opo ciona una es ela mucho menos oscilan e. En la
ale a mas co a obse amos como ambos ó ices se encuen an gene ando la oscilación de la es ela pe o con
meno in ensidad que en el caso sin ale a.
En el caso de la ale a mas la ga, igu a 4-26 obse amos como los o bellinos apenas es án en con ac o
gene ando una oscilación mínima de la es ela.
En la siguien e página mos amos el código MATLAB u ilizado pa a la esolución de es os p oblemas.
68
%=====================FLUJO ALREDEDOR DE UN PRISMA DE SECCIÓN
CUADRADA=================%
%================================PARÁMETROS=================================
===========%
close all;clea all;
nx=349;%in e alos de pa ición en x
ny=59;%in e alos de pa ición en y
n =20000;%in e alos de pa ición en
D=1; %ancho OBSTACULO
D =0.002; %inc emen o empo al
Dx=35*D/(nx); %inc emen o eje x
Dy=6*D/(ny); %inc emen o eje y
U=1; % elocidad ca idad
%===================================================================%
%===============NODOS DEL OBSTACULO=================================%
i1= ound(15/35*(nx+1));
i2= ound(16/35*(nx+1));
j1= ound(2.5/6*(ny+1));
j2= ound(3.5/6*(ny+1));
%===================================================================%
%===============NUMEROS ADIMENSIONALES==============================%
Re=40; %nume o de Reynolds
M=0.05; %nume o de Mach
%===================COEFICIENTES====================================%
a1=D /Dx;
a2=D /Dy;
a3=D /(Dx*M^2);
a4=D /(Dy*M^2);
a5=4*D /(3*Re*Dx^2);
a6=D /(Re*Dy^2);
a7=D /(Re*Dx^2);
a8=4*D /(3*Re*Dy^2);
a9=D /(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=ze os(nx+1,ny+1);% elocidad x en los pun os del mallado
P_ =ze os(nx+1,ny+1);% elocidad y en los pun os del mallado
P_ ho=ones(nx+1,ny+1);% densidad en los pun os del mallado
P_ hou=ze os(nx+1,ny+1);% densidad*u en los pun os del mallado
P_ ho =ze os(nx+1,ny+1);%densidad* en los pun os del mallado
P_ hou2=ze os(nx+1,ny+1);% densidad*u^2 en los pun os del mallado
P_ ho 2=ze os(nx+1,ny+1);%densidad* ^2 en los pun os del mallado
P_ hou =ze os(nx+1,ny+1);%densidad*u en los pun os del mallado
u=ze os(nx+1,ny+1);% elocidad x en los pun os del mallado
=ze os(nx+1,ny+1);% elocidad y en los pun os del mallado
ho=ones(nx+1,ny+1);% densidad en los pun os del mallado
hou=ze os(nx+1,ny+1);% densidad*u en los pun os del mallado
ho =ze os(nx+1,ny+1);%densidad* en los pun os del mallado
hou2= hou.*u(:,:);% densidad*u^2 en los pun os del mallado
ho 2=ze os(nx+1,ny+1);%densidad* ^2 en los pun os del mallado
hou =ze os(nx+1,ny+1);%densidad*u en los pun os del mallado
%========================BUCLE TEMPORAL===============================%
con =1;
k=0;
o q=1:n
69
k=k+1;
q
%=========================PREDICTOR==================================%
P_ ho(2:nx,2:ny)= ho(2:nx,2:ny)-a1*( hou(3:nx+1,2:ny)- hou(2:nx,2:ny))-
a2*( ho (2:nx,3:ny+1)- ho (2:nx,2:ny));
P_ hou(2:nx,2:ny)= hou(2:nx,2:ny)-a3*( ho(3:nx+1,2:ny)- ho(2:nx,2:ny))-
a1*( hou2(3:nx+1,2:ny)- hou2(2:nx,2:ny))-a2*( hou (2:nx,3:ny+1)-
hou (2:nx,2:ny))-a10*u(2:nx,2:ny)+a5*(u(3:nx+1,2:ny)+u(1:nx-
1,2:ny))+a6*(u(2:nx,3:ny+1)+u(2:nx,1:ny-1))+a9*( (3:nx+1,3:ny+1)+ (1:nx-
1,1:ny-1)- (1:nx-1,3:ny+1)- (3:nx+1,1:ny-1));
P_ ho (2:nx,2:ny)= ho (2:nx,2:ny)-a4*( ho(2:nx,3:ny+1)- ho(2:nx,2:ny))-
a1*(( hou (3:nx+1,2:ny))-( hou (2:nx,2:ny)))-a2*( ho 2(2:nx,3:ny+1)-
ho 2(2:nx,2:ny))-a11* (2:nx,2:ny)+a7*( (3:nx+1,2:ny)+ (1:nx-
1,2:ny))+a8*( (2:nx,3:ny+1)+ (2:nx,1:ny-1))+a9*(u(3:nx+1,3:ny+1)+u(1:nx-
1,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_ hou(2:nx,2:ny)./P_ ho(2:nx,2:ny);
P_ (2:nx,2:ny)=P_ ho (2:nx,2:ny)./P_ ho(2:nx,2:ny);
%=============================C.C====================================%
P_ ho(1,1:ny+1)= ho(1,1:ny+1)-(a1/2)*(- hou(3,1:ny+1)+4* hou(2,1:ny+1)-
3* hou(1,1:ny+1));%%%%x=0 densidad
P_ ho(nx+1,1:ny+1)= ho(nx+1,1:ny+1)+(a1/2)*(- hou(nx-
1,1:ny+1)+4* hou(nx,1:ny+1)-3* hou(nx+1,1:ny+1));%%%x=d densidad
P_ ho(2:nx,1)= ho(2:nx,1)-(a2/2)*(- ho (2:nx,3)+4* ho (2:nx,2)-
3* ho (2:nx,1))-(a1*U/2)*( ho(3:nx+1,1)- ho(1:nx-1,1));%%%%y=0 densidad
P_ ho(2:nx,ny+1)= ho(2:nx,ny+1)+a2/2*(- ho (2:nx,ny-1)+4* ho (2:nx,ny)-
3* ho (2:nx,ny+1))-(a1*U/2)*( ho(3:nx+1,ny+1)- ho(1:nx-1,ny+1));%%%y=d
densidad
P_ ho(i1,j1:j2)=1/3*(4* ho(i1-1,j1:j2)- ho(i1-2,j1:j2))+a12*(-5*u(i1-
1,j1:j2)+4*u(i1-2,j1:j2)-u(i1-3,j1:j2))-a13*(-( (i1-2,j1+1:j2+1)- (i1,j1-
1:j2-1))+4*( (i1-1,j1+1:j2+1)- (i1-1,j1-1:j2-1))-3*( (i1,j1+1:j2+1)- (i1,j1-
1:j2-1)));%%%%x=0 densidad cubo
P_ ho(i2,j1:j2)=1/3*(4* ho(i2+1,j1:j2)- ho(i2+2,j1:j2))-a12*(-
5*u(i2+1,j1:j2)+4*u(i2+2,j1:j2)-u(i2+3,j1:j2))-a13*(-( (i2+2,j1+1:j2+1)-
(i2,j1-1:j2-1))+4*( (i2+1,j1+1:j2+1)- (i2+1,j1-1:j2-1))-3*( (i2,j1+1:j2+1)-
(i2,j1-1:j2-1)));%%%x=d densidad cubo
P_ ho(i1:i2,j1)=1/3*(4* ho(i1:i2,j1-1)- ho(i1:i2,j1-2))+a12*(-
5* (i1:i2,j1-1)+4* (i1:i2,j1-2)- (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_ ho(i1:i2,j2)=1/3*(4* ho(i1:i2,j2+1)- ho(i1:i2,j2+2))-a12*(-
5* (i1:i2,j2+1)+4* (i1:i2,j2+2)- (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;%%%% elocidad u y=d
P_ (2:nx,ny+1)=0;%%%% elocidad y=d
P_u(2:nx,1)=U;%%%% elocidad u y=0
P_ (2:nx,1)=0;%%%% elocidad y=0
P_u(nx+1,:)=U;%%%% elocidad u x=d
P_ (nx+1,:)=0;%%%% elocidad x=d
P_u(1,:)=U;%%%% elocidad u x=0
P_ (1,:)=0;%%%% elocidad x=0
P_u(i1:i2,j1:j2)=0;
P_ (i1:i2,j1:j2)=0;
%===========================ACTUALIZACIÓN========================%
P_ hou=P_ ho.*P_u;
P_ ho =P_ ho.*P_ ;
P_ hou2=P_ hou.*P_u;
P_ hou =P_ hou.*P_ ;
P_ ho 2=P_ ho .*P_ ;
%==========================CORRECTOR==============================%
ho(2:nx,2:ny)=( ho(2:nx,2:ny)+P_ ho(2:nx,2:ny) -a1*(P_ hou(2:nx,2:ny)-
P_ hou(1:nx-1,2:ny))-a2*(P_ ho (2:nx,2:ny)-P_ ho (2:nx,1:ny-1)))*0.5;
hou(2:nx,2:ny)=0.5*( hou(2:nx,2:ny)+P_ hou(2:nx,2:ny)-
a3*(P_ ho(2:nx,2:ny)-P_ ho(1:nx-1,2:ny))-a1*(P_ hou2(2:nx,2:ny)-
P_ hou2(1:nx-1,2:ny))-a2*(P_ hou (2:nx,2:ny)-P_ hou (2:nx,1:ny-1))-
a10*P_u(2:nx,2:ny)+a5*(P_u(3:nx+1,2:ny)+P_u(1:nx-
1,2:ny))+a6*(P_u(2:nx,3:ny+1)+P_u(2:nx,1:ny-
1))+a9*(P_ (3:nx+1,3:ny+1)+P_ (1:nx-1,1:ny-1)-P_ (1:nx-1,3:ny+1)-
P_ (3:nx+1,1:ny-1)));
ho (2:nx,2:ny)= 0.5*( ho (2:nx,2:ny)+P_ ho (2:nx,2:ny)-
a4*(P_ ho(2:nx,2:ny)-P_ ho(2:nx,1:ny-1))-a1*(P_ hou (2:nx,2:ny)-
P_ hou (1:nx-1,2:ny))-a2*(P_ ho 2(2:nx,2:ny)-P_ ho 2(2:nx,1:ny-1))-
a11*P_ (2:nx,2:ny)+a7*(P_ (3:nx+1,2:ny)+P_ (1:nx-
1,2:ny))+a8*(P_ (2:nx,3:ny+1)+P_ (2:nx,1:ny-
1))+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= hou./ ho;
= ho ./ ho;
%=============================C.C====================================%
ho(1,2:ny)=(P_ ho(1,2:ny)+ ho(1,2:ny)-a1/2*(-
P_ hou(3,2:ny)+4*P_ hou(2,2:ny)-3*P_ hou(1,2:ny)))*0.5;%%%%x=0
ho(nx+1,2:ny)=(P_ ho(nx+1,2:ny)+ ho(1,2:ny)+a1/2*(-P_ hou(nx-
1,2:ny)+4*P_ hou(nx,2:ny)-3*P_ hou(nx+1,2:ny)))*0.5;%%%x=d
ho(2:nx,1)=(P_ ho(2:nx,1)+ ho(2:nx,1)-a2/2*(-
P_ ho (2:nx,3)+4*P_ ho (2:nx,2)-3*P_ ho (2:nx,1))-a1*U/2*(P_ ho(3:nx+1,1)-
P_ ho(1:nx-1,1)))*0.5;%%%%y=0
ho(2:nx,1+ny)=(P_ ho(2:nx,1+ny)+ ho(2:nx,1+ny)+a2/2*(-P_ ho (2:nx,ny-
1)+4*P_ ho (2:nx,ny)-3*P_ ho (2:nx,ny+1))-a1*U/2*(P_ ho(3:nx+1,ny+1)-
P_ ho(1:nx-1,ny+1)))*0.5;%%%y=d
71
ho(i1,j1:j2)=( ho(i1,j1:j2)+1/3*(4*P_ ho(i1-1,j1:j2)-P_ ho(i1-
2,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_ (i1-2,j1+1:j2+1)-P_ (i1,j1-1:j2-1))+4*(P_ (i1-1,j1+1:j2+1)-P_ (i1-1,j1-
1:j2-1))-3*(P_ (i1,j1+1:j2+1)-P_ (i1,j1-1:j2-1))))*0.5;%%%%x=0 densidad cubo
ho(i2,j1:j2)=( ho(i2,j1:j2)+1/3*(4*P_ ho(i2+1,j1:j2)-
P_ ho(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_ (i2+2,j1+1:j2+1)-P_ (i2,j1-1:j2-
1))+4*(P_ (i2+1,j1+1:j2+1)-P_ (i2+1,j1-1:j2-1))-3*(P_ (i2,j1+1:j2+1)-
P_ (i2,j1-1:j2-1))))*0.5;%%%x=d densidad cubo
ho(i1:i2,j1)=( ho(i1:i2,j1)+1/3*(4*P_ ho(i1:i2,j1-1)-P_ ho(i1:i2,j1-
2))+a12*(-5*P_ (i1:i2,j1-1)+4*P_ (i1:i2,j1-2)-P_ (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(i1-
1: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
ho(i1:i2,j2)=( ho(i1:i2,j2)+1/3*(4*P_ ho(i1:i2,j2+1)-
P_ ho(i1:i2,j2+2))-a12*(-5*P_ (i1:i2,j2+1)+4*P_ (i1:i2,j2+2)-
P_ (i1:i2,j2+3))-a13*(-(P_u(i1+1:i2+1,j2+2)-P_u(i1-1:i2-
1,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;%%%% elocidad u y=d
(2:nx,ny+1)=0;%%%% elocidad y=d
u(2:nx,1)=U;%%%% elocidad u y=0
(2:nx,1)=0;%%%% elocidad y=0
u(nx+1,:)=U;%%%% elocidad u x=d
(nx+1,:)=0;%%%% elocidad x=d
u(1,:)=U;%%%% elocidad u x=0
(1,:)=0;%%%% elocidad x=0
u(i1:i2,j1:j2)=0;
(i1:i2,j1:j2)=0;
%===========================ACTUALIZACIÓN========================%
hou= ho.*u;
ho = ho.* ;
hou2= hou.*u;
hou = hou.* ;
ho 2= ho .* ;
%======================REPRESENTACIÓN GRÁFICA====================%
i k==100
pel(:,:,con )=[u ];
k=0;
con =con +1;
end
%======================CÁLCULO DE COEFICIENTES===================%
CL(q)=2*sum( ho(i1:i2,j1)- ho(i1:i2,j2))*Dx/M^2;
CD(q)=2*sum( ho(i1,j1:j2)- ho(i2,j1:j2))*Dy/M^2;
72
end
igu e
T=sq ( '.^2+u'.^2);
qui e (linspace(0,35,nx+1),linspace(0,6,ny+1),u'./T, './T)
=D *n
73
5 APLICACIÓN: FLUJO
BIDIMENSIONAL ALREDEDOR DE
UN CUADRADO CON
TEMPERATURA
El esquema de es a aplicación es el mismo que la an e io , un cuad ado desplazándose hacia la izquie da en
un canal, pe o en es e caso la empe a u a del cuad ado se á mayo a la empe a u a ambien e pe mi iéndonos
asi es udia las capas limi e Té micas que pudiesen gene a se.
Figu a( 5-1 )
Las ecuaciones en es a aplicación lógicamen e deben de se complemen adas con la ecuación de la ene gía y
condicionesde con o no en empe a u a. A con inuación se mues an las ecuaciones comple as que
expusimos en el apa ado 3 de las que ealiza emos algunas simpli icaciones:
𝜕(𝑈)
𝜕𝑡 +𝜕(𝐹)
𝜕𝑥 +𝜕(𝐺)
𝜕𝑦 +𝑆=0
( 5-1)
𝑈=[𝜌
𝜌𝑢
𝜌𝑣
𝐸]
( 5-2)
𝐹=
[
𝜌𝑢
𝜌𝑢2+𝑝−𝜏𝑥𝑥
𝜌𝑢𝑣−𝜏𝑥𝑦
(𝐸+𝑝)𝑢−𝑢𝜏𝑥𝑥−𝑣𝜏𝑥𝑦+𝑞𝑥
]
( 5-3)
80
Donde se ap ecia con más cla idad si cabe la dominancia de los mecanismos de ans e encia de calo po
con ección en el caso de 𝑃𝑟=7 mien as que en el caso de 𝑃𝑟=0.7 ambos mecanismos son impo an es.
Si ealizamos un análisis de ó denes de magni ud podemos es ima los espeso es de la capa lími e é mica.
𝛿𝑡
𝐷=1
√𝑅𝑒𝑃𝑟
( 5-34)
Teniendo en cuen a que ambos caos se han ob enido pa a un mismo núme o de Reynolds podemos ob ene
la siguien e elación:
𝛿𝑡0.7
𝛿𝑡7=√7
0.7=√10~3.16
( 5-35)
Es o se puede co obo a obse ando los pe iles de empe a u a ep esen ados an e io men e:
𝛿𝑡0.7
𝛿𝑡7=2.5
0.9375=2.668
( 5-36)
Vemos como ambos núme os son del mismo o den de magni ud. No obs an e, es con enien e comen a que
en el caso de 𝑃𝑟=0.7 los e ec os pe u bado es de la pa ed del conduc o apa ecen an e de que se pueda
desa olla comple amen e el pe il de empe a u as.
En la página siguien e mos amos el código MATLAB u ilizado pa a la esolución de es os p oblemas.

81
%=======FLUJO ALREDEDOR DE UN CUADRADO CON TEMPERATURA=========%
%LAS VARIABLES P_ SE CORRESPONDEN AL TERMINO PREDICTOR
%==============================================================%
close all;clea all;
%========================PARÁMETROS============================%
nx=319;%in e alos de pa ición en x
ny=179;%in e alos de pa ición en y
n =35000;%in e alos de pa ición en
D=1; %ancho ca idad
D =0.003; %inc emen o empo al
Dx=35*D/nx; %inc emen o eje x
Dy=6*D/ny; %inc emen o eje y
%=======================NODOS CUADRADO==========================%
i1= ound(15/35*(nx+1));
i2= ound(16/35*(nx+1));
j1= ound(2.5/6*(ny+1));
j2= ound(3.5/6*(ny+1));
%========================NÚMEROS ADIMENSIONALES=================%
M=0.15; %núme o de Mach
Re=100; %núme o de Reynolds
gamma=1.4;
d =0.35; %inc emen o empe a u a
%=========================TEMPERATURAS Y
VELOCIDAD==========================%
T2=1+d ;% empe a u a cuad ado
T1=1; % empe a u a luido
U=M; % elocidad ca idad
%==========================COEFICIENTES==================================%
a1=D /Dx;
a2=D /Dy;
a3=D /F ^2;
a4=-D /P /Re/(gamma-1)/Dx^2;
a52=-D /P /Re/(gamma-1)/Dy^2;
a5=4*D /(3*Re*Dx^2);
a6=D /(Re*Dy^2);
a7=D /(Re*Dx^2);
a8=4*D /(3*Re*Dy^2);
a9=D /(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]=meshg id(x,y) ; %mallado
%==================CONDICIONES INICIALES=============%
P_u=ones(nx+1,ny+1)*U;% elocidad x en los pun os del mallado
P_ =ze os(nx+1,ny+1);% elocidad y en los pun os del mallado
P_u(i1:i2,j1:j2)=0; % elocidad x en el cuad ado
P_ (i1:i2,j1:j2)=0; % elocidad y en el cuad ado
P_T=ones(nx+1,ny+1)*T1; %campo de empe a u a en cada pun o del mallado
equilib io
P_ ho=ones(nx+1,ny+1);% densidad en los pun os del mallado
P_ hou=P_ ho.*P_u;% densidad*u en los pun os del mallado
P_ ho =P_ ho.*P_ ;%densidad* en los pun os del mallado
P_ hou2=P_ hou.*P_u;% densidad*u^2 en los pun os del mallado
P_ ho 2=P_ ho .*P_ ;%densidad* ^2 en los pun os del mallado
82
P_ hou =P_ hou.*P_ ;%densidad*u en los pun os del mallado
P_E=(P_T/gamma/(gamma-1)-0.5*(P_ .^2+P_u.^2)).*P_ ho;%ene gía o al en los
pun os del mallado
P_Eu=P_E.*P_u;%ene gía*u o al en los pun os del mallado
P_E =P_E.*P_ ;%ene gía* o al en los pun os del mallado
u=ones(nx+1,ny+1)*U;% elocidad x en los pun os del mallado
=ze os(nx+1,ny+1);% elocidad y en los pun os del mallado
u(i1:i2,j1:j2)=0; % elocidad x en el cuad ado
(i1:i2,j1:j2)=0;% elocidad en el cuad ado
T=ones(nx+1,ny+1)*T1; %campo de empe a u a en cada pun o del mallado
equilib io
ho=ones(nx+1,ny+1);% densidad en los pun os del mallado
hou= ho.*u;% densidad*u en los pun os del mallado
ho = ho.* ;%densidad* en los pun os del mallado
hou2= hou.*u;% densidad*u^2 en los pun os del mallado
ho 2= ho .* ;%densidad* ^2 en los pun os del mallado
hou = hou.* ;%densidad*u en los pun os del mallado
E=(T/gamma/(gamma-1)-0.5*( .^2+u.^2)).* ho;%ene gía o al en los pun os del
mallado
Eu=E.*u;%ene gía*u o al en los pun os del mallado
E =E.* ;%ene gía* o al en los pun os del mallado
%==================BUCLE TEMPORAL===========================%
con =1;
k=1;
o q=1:n
u(:,:)= hou(:,:)./ ho(:,:);
(:,:)= ho (:,:)./ ho(:,:);
k=1+k;
%===========================PREDICTOR=============================%
pdgamma=(gamma-1)*(E-0.5* ho.*(u.^2+ .^2));%VARIABLE AUXILIAR PRESIÓN
P_ ho(2:nx,2:ny)= ho(2:nx,2:ny)-a1*( hou(3:nx+1,2:ny)- hou(2:nx,2:ny))-
a2*( ho (2:nx,3:ny+1)- ho (2:nx,2:ny));
P_ hou(2:nx,2:ny)= hou(2:nx,2:ny)-a1*(gamma-1)*(E(3:nx+1,2:ny)-
E(2:nx,2:ny))+a1*0.5*(gamma-1)*(( ho 2(3:nx+1,2:ny)+ hou2(3:nx+1,2:ny))-
( ho 2(2:nx,2:ny)+ hou2(2:nx,2:ny)))-a1*( hou2(3:nx+1,2:ny)-
hou2(2:nx,2:ny))-a2*( hou (2:nx,3:ny+1)- hou (2:nx,2:ny))-
a10*u(2:nx,2:ny)+a5*(u(3:nx+1,2:ny)+u(1:nx-
1,2:ny))+a6*(u(2:nx,3:ny+1)+u(2:nx,1:ny-1))+a9*( (3:nx+1,3:ny+1)+ (1:nx-
1,1:ny-1)- (1:nx-1,3:ny+1)- (3:nx+1,1:ny-1));
P_ ho (2:nx,2:ny)= ho (2:nx,2:ny)-a2*(gamma-1)*(E(2:nx,3:ny+1)-
E(2:nx,2:ny))+a1*0.5*(gamma-1)*(( ho 2(2:nx,3:ny+1)+ hou2(2:nx,3:ny+1))-
( ho 2(2:nx,2:ny)+ hou2(2:nx,2:ny)))-a1*(( hou (3:nx+1,2:ny))-
( hou (2:nx,2:ny)))-a2*( ho 2(2:nx,3:ny+1)- ho 2(2:nx,2:ny))-
a11* (2:nx,2:ny)+a7*( (3:nx+1,2:ny)+ (1:nx-
1,2:ny))+a8*( (2:nx,3:ny+1)+ (2:nx,1:ny-1))+a9*(u(3:nx+1,3:ny+1)+u(1:nx-
1,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*( (2:nx,3:ny+1).*(E(2:nx,3:ny+1)+pdgamma(2:nx,3:ny+1))-
(2:nx,2:ny).*(E(2:nx,2:ny)+pdgamma(2:nx,2:ny)));
83
%=============================DECODIFICACIÓN VARIABLES==================%
P_u=P_ hou./P_ ho;
P_ =P_ ho ./P_ ho;
Vmod2(2:nx,2:ny)=P_ (2:nx,2:ny).^2+P_u(2:nx,2:ny).^2;% a iable auxilia
módulo elocidad al cuad ado
P_T(2:nx,2:ny)=(P_E(2:nx,2:ny)./P_ ho(2:nx,2:ny)-
0.5*Vmod2(2:nx,2:ny))*gamma*(gamma-1);
%========================CONDICIONES DE CONTORNO PREDICTOR==============%
P_ ho(1,1:ny+1)= ho(1,1:ny+1)-(a1/2)*(- hou(3,1:ny+1)+4* hou(2,1:ny+1)-
3* hou(1,1:ny+1));%%%%x=0 densidad
P_ ho(nx+1,1:ny+1)= ho(nx+1,1:ny+1)+(a1/2)*(- hou(nx-
1,1:ny+1)+4* hou(nx,1:ny+1)-3* hou(nx+1,1:ny+1));%%%x=d densidad
P_ ho(2:nx,1)= ho(2:nx,1)-(a2/2)*(- ho (2:nx,3)+4* ho (2:nx,2)-
3* ho (2:nx,1))-(a1*U/2)*( ho(3:nx+1,1)- ho(1:nx-1,1));%%%%y=0 densidad
P_ ho(2:nx,ny+1)= ho(2:nx,ny+1)+a2/2*(- ho (2:nx,ny-1)+4* ho (2:nx,ny)-
3* ho (2:nx,ny+1))-(a1*U/2)*( ho(3:nx+1,ny+1)- ho(1:nx-1,ny+1));%%%y=d
densidad
P_ ho(i1,j1:j2)= ho(i1,j1:j2)+(a1/2)*(- hou(i1-2,j1:j2)+4* hou(i1-
1,j1:j2)-3* hou(i1,j1:j2));%%%%x=0 densidad cubo
P_ ho(i2,j1:j2)= ho(i2,j1:j2)-(a1/2)*(-
hou(i2+2,j1:j2)+4* hou(i2+1,j1:j2)-3* hou(i2,j1:j2));%%%x=d densidad cubo
P_ ho(i1:i2,j1)= ho(i1:i2,j1)+a2/2*(- ho (i1:i2,j1-2)+4* ho (i1:i2,j1-1)-
3* ho (i1:i2,j1));%%%%y=0 densidad cubo
P_ ho(i1:i2,j2)= ho(i1:i2,j2)-(a2/2)*(-
ho (i1:i2,j2+2)+4* ho (i1:i2,j2+1)-3* ho (i1:i2,j2));%%%y=d densidad cubo
P_u(2:nx,ny+1)=U;%%%% elocidad u y=d
P_ (2:nx,ny+1)=0;%%%% elocidad y=d
P_u(2:nx,1)=U;%%%% elocidad u y=0
P_ (2:nx,1)=0;%%%% elocidad y=0
P_u(nx+1,:)=U;%%%% elocidad u x=d
P_ (nx+1,:)=0;%%%% elocidad x=d
P_u(1,:)=U;%%%% elocidad u x=0
P_ (1,:)=0;%%%% elocidad x=0
P_u(i1:i2,j2)=0;%%%% elocidad u y=d cubo
P_ (i1:i2,j2)=0;%%%% elocidad y=d cubo
P_u(i1:i2,j1)=0;%%%% elocidad u y=0 cubo
P_ (i1:i2,j1)=0;%%%% elocidad y=0 cubo
P_u(i2,j1:j2)=0;%%%% elocidad u x=d cubo
P_ (i2,j1:j2)=0;%%%% elocidad x=d cubo
P_u(i1,j1:j2)=0;%%%% elocidad u x=0 cubo
P_ (i1,j1:j2)=0;%%%% elocidad x=0 cubo
P_u(i1:i2,j1:j2)=0;
P_ (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_ ho(:,1)/gamma/(gamma-1);
P_E(:,ny+1)=P_T(:,ny+1).*P_ ho(:,ny+1)/gamma/(gamma-1);
P_E(1,2:ny)=P_T(1,2:ny).*P_ ho(1,2:ny)/gamma/(gamma-1);
P_E(nx+1,2:ny)=P_T(nx+1,2:ny).*P_ ho(nx+1,2:ny)/gamma/(gamma-1);
P_ hou(2:nx,ny+1)=P_ ho(2:nx,1+ny).* P_u(2:nx,ny+1);%%%%%%%densidad*u en
y=d
P_ hou(2:nx,1)=P_ ho(2:nx,1).* P_u(2:nx,1);%%%%%%%densidad*u en y=0
P_ hou(1,:)=P_ ho(1,:).* P_u(1,:);%%%%%%%densidad*u en x=0
P_ hou(nx+1,:)=P_ ho(nx+1,:).* P_u(nx+1,:);%%%%%%%densidad*u en x=d
P_ ho (2:nx,ny+1)=P_ ho(2:nx,1+ny).* P_ (2:nx,ny+1);%%%%%%%densidad* en
y=d
P_ ho (2:nx,1)=P_ ho(2:nx,1).* P_ (2:nx,1);%%%%%%%densidad* en y=0
P_ ho (1,:)=P_ ho(1,:).* P_ (1,:);%%%%%%%densidad* en x=0
P_ ho (nx+1,:)=P_ ho(nx+1,:).* P_ (nx+1,:);%%%%%%%densidad* en x=d
%=======================ACTUALIZACIÓN DE VARIABLES=============%
Vmod2(2:nx,2:ny)=P_ (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_ ho(2:nx,2:ny)-
0.5*Vmod2(2:nx,2:ny))*gamma*(gamma-1);
P_ hou2=P_ hou.*P_u;
P_ hou =P_ hou.*P_ ;
P_ ho 2=P_ ho .*P_ ;
P_Eu=P_E.*P_u;
P_E =P_E.*P_ ;
%============================CORRECTOR========================%
pdgamma=0.5*(pdgamma+(gamma-1).*(P_E-0.5.*P_ ho.*(P_u.^2+P_ .^2)));
ho(2:nx,2:ny)=0.5*(P_ ho(2:nx,2:ny)+ ho(2:nx,2:ny)-a1*( hou(2:nx,2:ny)-
hou(1:nx-1,2:ny))-a2*( ho (2:nx,2:ny)- ho (2:nx,1:ny-1)));
hou(2:nx,2:ny)=0.5*(P_ hou(2:nx,2:ny)+ hou(2:nx,2:ny)-a1*(gamma-
1)*(P_E(2:nx,2:ny)-P_E(1:nx-1,2:ny))+a1*0.5*(gamma-
1)*(P_ ho 2(2:nx,2:ny)+P_ hou2(2:nx,2:ny)-(P_ ho 2(1:nx-
1,2:ny)+P_ hou2(1:nx-1,2:ny)))-a1*( hou2(2:nx,2:ny)- hou2(1:nx-1,2:ny))-
a2*(P_ hou (2:nx,2:ny)-P_ hou (2:nx,1:ny-1))-
a10*P_u(2:nx,2:ny)+a5*(P_u(3:nx+1,2:ny)+P_u(1:nx-
1,2:ny))+a6*(P_u(2:nx,3:ny+1)+P_u(2:nx,1:ny-
1))+a9*(P_ (3:nx+1,3:ny+1)+P_ (1:nx-1,1:ny-1)-P_ (1:nx-1,3:ny+1)-
P_ (3:nx+1,1:ny-1)));
ho (2:nx,2:ny)=0.5*(P_ ho (2:nx,2:ny)+ ho (2:nx,2:ny)-a2*(gamma-
1)*(P_E(2:nx,2:ny)-P_E(2:nx,1:ny-1))+a2*0.5*(gamma-
1)*(P_ ho 2(2:nx,2:ny)+P_ hou2(2:nx,2:ny)-(P_ ho 2(2:nx,1:ny-
1)+P_ hou2(2:nx,1:ny-1)))-a1*((P_ hou (2:nx,2:ny))-(P_ hou (1:nx-1,2:ny)))-
a2*( ho 2(2:nx,2:ny)- ho 2(2:nx,1:ny-1))-
a11* (2:nx,2:ny)+a7*( (3:nx+1,2:ny)+ (1:nx-
1,2:ny))+a8*( (2:nx,3:ny+1)+ (2:nx,1:ny-1))+a9*(u(3:nx+1,3:ny+1)+u(1:nx-
1,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:nx-
1,2:ny)+pdgamma(1:nx-1,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*( (2:nx,2:ny).*(E(2:nx,2:ny)+pdgamma(2:nx,2:ny))- (2:nx,1:ny-
1).*(E(2:nx,1:ny-1)+pdgamma(2:nx,1:ny-1))));
%================ACTUALIZACIÓN DE VARIABLES=====================%
u= hou./ ho;
= ho ./ ho;
Vmod2(2:nx,2:ny)= (2:nx,2:ny).^2+u(2:nx,2:ny).^2;
85
T(2:nx,2:ny)=(E(2:nx,2:ny)./ ho(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 pa ed calien e
T(:,ny+1)=(4*T(:,ny)-T(:,ny-1))/3;%%%%y=d pa ed ia
T(1,2:ny)=T1 ; % T(nx+1,2:ny)=T(nx,2:ny)+(T(nx,2:ny)-T(nx-
1,2:ny));
T(i1:i2,j2)=T2;%%%% elocidad u y=d cubo
T(i1:i2,j1)=T2;%%%% elocidad u y=0 cubo
T(i2,j1:j2)=T2;%%%% elocidad u x=d cubo
T(i1,j1:j2)=T2;%%%% elocidad u x=0 cubo
T(i1:i2,j1:j2)=T2;
ho(1,2:ny)=(P_ ho(1,2:ny)+ ho(1,2:ny)-a1/2*(-
P_ hou(3,2:ny)+4*P_ hou(2,2:ny)-3*P_ hou(1,2:ny)))*0.5;%%%%x=0
ho(nx+1,2:ny)=(P_ ho(nx+1,2:ny)+ ho(1,2:ny)+a1/2*(-P_ hou(nx-
1,2:ny)+4*P_ hou(nx,2:ny)-3*P_ hou(nx+1,2:ny)))*0.5;%%%x=d
ho(2:nx,1)=(P_ ho(2:nx,1)+ ho(2:nx,1)-a2/2*(-
P_ ho (2:nx,3)+4*P_ ho (2:nx,2)-3*P_ ho (2:nx,1))-a1*U/2*(P_ ho(3:nx+1,1)-
P_ ho(1:nx-1,1)))*0.5;%%%%y=0
ho(2:nx,1+ny)=(P_ ho(2:nx,1+ny)+ ho(2:nx,1+ny)+a2/2*(-P_ ho (2:nx,ny-
1)+4*P_ ho (2:nx,ny)-3*P_ ho (2:nx,ny+1))-a1*U/2*(P_ ho(3:nx+1,ny+1)-
P_ ho(1:nx-1,ny+1)))*0.5;%%%y=d
ho(i1,j1:j2)=( ho(i1,j1:j2)+P_ ho(i1,j1:j2)+(a1/2)*(-P_ hou(i1-
2,j1:j2)+4*P_ hou(i1-1,j1:j2)-3*P_ hou(i1,j1:j2)))*0.5;%%%%x=0 densidad cubo
ho(i2,j1:j2)=( ho(i2,j1:j2)+P_ ho(i2,j1:j2)-(a1/2)*(-
P_ hou(i2+2,j1:j2)+4*P_ hou(i2+1,j1:j2)-3*P_ hou(i2,j1:j2)))*0.5;%%%x=d
densidad cubo
ho(i1:i2,j1)=( ho(i1:i2,j1)+P_ ho(i1:i2,j1)+a2/2*(-P_ ho (i1:i2,j1-
2)+4*P_ ho (i1:i2,j1-1)-3*P_ ho (i1:i2,j1)))*0.5;%%%%y=0 densidad cubo
ho(i1:i2,j2)=( ho(i1:i2,j2)+P_ ho(i1:i2,j2)-(a2/2)*(-
P_ ho (i1:i2,j2+2)+4*P_ ho (i1:i2,j2+1)-3*P_ ho (i1:i2,j2)))*0.5;%%%y=d
densidad cubo
u(2:nx,ny+1)=U;%%%% elocidad u y=d
(2:nx,ny+1)=0;%%%% elocidad y=d
u(2:nx,1)=U;%%%% elocidad u y=0
(2:nx,1)=0;%%%% elocidad y=0
u(nx+1,:)=U;%%%% elocidad u x=d
(nx+1,:)=0;%%%% elocidad x=d
u(1,:)=U;%%%% elocidad u x=0
(1,:)=0;%%%% elocidad x=0
u(i1:i2,j2)=0;%%%% elocidad u y=d cubo
(i1:i2,j2)=0;%%%% elocidad y=d cubo
u(i1:i2,j1)=0;%%%% elocidad u y=0 cubo
(i1:i2,j1)=0;%%%% elocidad y=0 cubo
u(i2,j1:j2)=0;%%%% elocidad u x=d cubo
(i2,j1:j2)=0;%%%% elocidad x=d cubo
u(i1,j1:j2)=0;%%%% elocidad u x=0 cubo
(i1,j1:j2)=0;%%%% elocidad x=0 cubo
u(i1:i2,j1:j2)=0;
(i1:i2,j1:j2)=0;

86
E(:,1)=T(:,1).* ho(:,1)/gamma/(gamma-1);
E(:,ny+1)=T(:,ny+1).* ho(:,ny+1)/gamma/(gamma-1);
E(1,2:ny)=T(1,2:ny).* ho(1,2:ny)/gamma/(gamma-1);
E(nx+1,2:ny)=T(nx+1,2:ny).* ho(nx+1,2:ny)/gamma/(gamma-1);
hou(2:nx,ny+1)= ho(2:nx,1+ny).* u(2:nx,ny+1);%%%%%%%densidad*u en y=d
hou(2:nx,1)= ho(2:nx,1).*u(2:nx,1);%%%%%%%densidad*u en y=0
hou(1,:)= ho(1,:).* u(1,:);%%%%%%%densidad*u en x=0
hou(nx+1,:)= ho(nx+1,:).*u(nx+1,:);%%%%%%%densidad*u en x=d
ho (2:nx,ny+1)= ho(2:nx,1+ny).* (2:nx,ny+1);%%%%%%%densidad* en y=d
ho (2:nx,1)= ho(2:nx,1).* (2:nx,1);%%%%%%%densidad* en y=0
ho (1,:)= ho(1,:).* (1,:);%%%%%%%densidad* en x=0
ho (nx+1,:)= ho(nx+1,:).* (nx+1,:);%%%%%%%densidad* en x=d
%===================ACTUALIZACIÓN DE VARIABLES=====================%
Eu=E.*u;
E =E.* ;
hou= ho.*u;% densidad*u en los pun os del mallado
ho = ho.* ;%densidad* en los pun os del mallado
hou2= hou.*u;
hou = hou.* ;
ho 2= ho .* ;
q;
%=============REPRESENTACIÓN GRÁFICA==========================%
i k==50
disp(max(max(abs( ))));
% plo (H)
% pause
q
% igu e
% con ou (T')
% Tmod=sq ( '.^2+u'.^2);
% igu e
% con ou (Tmod)
% igu e
% qui e (linspace(0,2,nx+1),linspace(0,1,ny+1),u'./Tmod, './Tmod)
k=1
pel(:,:,con )=T;
con =con +1;
end
H(q)=max(max(abs(u)));
end
qui e (linspace(0,2,nx+1),linspace(0,1,ny+1),u'./sq ( '.^2+u'.^2), './sq (
'.^2+u'.^2))
87
6 APLICACIÓN: CAVIDAD
RECTANGULAR CON DIFERENCIA
DE TEMPERATURA EN LAS TAPAS
El obje i o de es a aplicación se á ep oduci los pa ones de con ección na u al o iginados po la
ines abilidad que gene an las ue zas de lo abilidad en un luido en eposo. Si enemos un luido en eposo y
se in oduce una pequeña pe u bación que p o oque un lige o desplazamien o e ical de la pa ícula, la
densidad media de la pa ícula se á meno que la de su en o no y end á una endencia aún más ascenden e.
De o ma opues a si el desplazamien o de la pa ícula es nega i o su densidad media se á mayo que la de su
en o no y es o p o oca á que su endencia sea aún más descenden e. Además de las ue zas de lo abilidad
apa ecen o as ue zas que ienden a amo igua el e ec o deses abilizado de la lo abilidad es as son los
e ec os de iscosidad y de conducción de calo . Muchos es udios de es e p oblema han sido ealizados
median e el uso de la ap oximación de Boussinesq y en el caso incomp esible, en es e manual se a a el
es udio con la esolución de las ecuaciones de Na ie -S okes comple as u ilizando un esquema de
MacCo mack, u ilizado po K.V.Pa che sky [19] pa a ealiza simulaciones a g an escala de la con ección
sola , po lo que se pod ía conside a como desa ollo u u o la aplicación de es as ecuaciones pa a la
esolución de p oblemas de ese ipo.
A con inuación, amos a ealiza la o mulación de la aplicación que aquí se a a esol e .
Figu a( 6-1 )
Se a a esol e el p oblema de con ección de Rayleigh-Bena d con la pa ed in e io a mayo empe a u a,
que como e emos es condición necesa ia pa a que se gene e la ines abilidad. Las pa edes la e ales amos a
conside a las pe ec amen e adiabá icas.
Comencemos con las ecuaciones que amos a esol e . Pa i emos de las ecuaciones de Na ie -S okes
comple as e i emos ealizando algunas simpli icaciones.
 Se añadi án ue zas másicas a las ecuaciones.
88
La elocidad ca ac e ís ica se puede es ima de iguala en la ecuación de can idad de mo imien o el
o den de magni ud de los é minos con ec i os y los de lo abilidad como:
𝑣𝑜=√Δ𝑇𝛽𝑔𝐻
( 6-1)
Además, el cuad ado del núme o de F oude es al que:
𝐹𝑟2=𝑔𝐻Δ𝑇𝛽
𝑔𝐻 ~𝑂(1)
( 6-2)
No podemos desp ecia las ue zas másicas lógicamen e ya que son las enca gadas de ealiza el
mo imien o del luido.
 Los é minos de disipación iscosa compa ados con los de con ección é mica se pueden exp esa
como:
𝜏:𝛻𝑣
𝜌𝑐𝑝𝑣 𝛻𝑇~𝜇𝑣02
𝐻2
𝜌𝑐𝑝𝑣0Δ𝑇
𝐻=𝜈
𝑐𝑝√𝑔𝐻Δ𝑇𝛽
𝐻Δ𝑇 ≪1
( 6-3)
Y po an o pod án se desp eciados los e ec os de disipación iscosa espec o a los de con ección
é mica.
𝜕(𝑈)
𝜕𝑡 +𝜕(𝐹)
𝜕𝑥 +𝜕(𝐺)
𝜕𝑦 +𝑆=0
( 6-4)
𝑈=[𝜌
𝜌𝑢
𝜌𝑣
𝐸]=[𝑈1
𝑈2
𝑈3
𝑈4]
( 6-5)
𝐹=[𝜌𝑢
𝜌𝑢2+𝑝−𝜏𝑥𝑥
𝜌𝑢𝑣−𝜏𝑥𝑦
(𝐸+𝑝)𝑢+𝑞𝑥]
( 6-6)
𝐺=
[
𝜌𝑣
𝜌𝑢𝑣−𝜏𝑥𝑦
𝜌𝑣2+𝑝−𝜏𝑦𝑦
(𝐸+𝑝)𝑣+𝑞𝑦
]
( 6-7)
𝑆=[00
𝜌𝑔
𝜌𝑣𝑔]
( 6-8)
Pa a e i a alo es muy al os de p esión y empe a u a y po an o educi el cálculo compu acional, amos a
ope a únicamen e con las di e encias de empe a u a y p esión espec o a las a iables en equilib io
hid os á ico:
𝑝=𝑝𝑒𝑞+𝑝+ , 𝑇=𝑇𝑒𝑞+𝑇+, 𝐸=𝐸𝑒𝑞+𝐸+
( 6-9)
89
Se in oducen las siguien es a iables adimensionales:
𝑥∗=𝑥𝐻, 𝑦∗=𝑦
𝐻,𝑣0=√ΔTβ𝑔𝐻=√ΔT∗𝑔𝐻 , 𝑢∗=𝑢
𝑣0 , 𝑣∗=𝑣
𝑣0, 𝑝∗
=𝑝+
𝜌0 02, 𝑝𝑒𝑞
∗=𝑝𝑒𝑞
𝜌0𝑅𝑔𝑇0, 𝜌∗=𝜌
𝜌0 , 𝑇∗=𝑇
𝑇0,
𝐸𝑒𝑞
∗=𝐸𝑒𝑞
𝜌0𝑅𝑔𝑇0, 𝐸∗=𝐸+
𝜌0 02
( 6-10)
Se de inen los núme os adimenionales:
𝑅𝑒=𝜌0𝑣0𝐻
𝜇 , 𝑃𝑟=𝜇𝑐𝑝
𝑘=𝜈𝛼 , 𝐹𝑟2=ΔTβ=Δ𝑇
𝑇0=Δ𝑇∗,𝑀=𝑣0
𝑎0
( 6-11)
Que ep esen an los Núme os de Reynolds, P and l, F oude y Mach espec i amen e. Siendo 𝑣0 una
elocidad ca ac e ís ica, 𝜇 la iscosidad del luido, 𝜌0 𝑦 𝑝0 la densidad y p esión de e e encia, 𝜅 la
conduc i idad é mica del luido, 𝑔 la acele ación g a i a o ia, H una dimensión ca ac e ís ica del p oblema
y 𝑎0 la elocidad del sonido.
Con es as nue as a iables adimensionales, de las que de aho a en adelan e se sup imi á el as e isco, el
sis ema de ecuaciones esul an e se á:
𝑑(𝑈)
𝑑𝑡 +𝑑(𝐹)
𝑑𝑥 +𝑑(𝐺)
𝑑𝑦 +𝑆=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
Figu a( 6-3 )
Figu a( 6-4 )

97
𝑅𝑎=4717
Figu a( 6-5 )
Figu a( 6-6 )
98
Obse amos como en el p ime caso, con el núme o de Rayleigh in e io al c í ico, el campo de empe a u as
se mues a sin al e a compa ado con la condición inicial. En el segundo caso el núme o de Rayleigh ha
supe ado el umb al a pa i del cual los e ec os de lo abilidad supe an a los iscosos y a los de conducción
de calo y el luido comienza a ascende po el cen o y a descende po las pa edes o mando es e pa ón
celula an ep esen a i o. Es e pa ón de lujo p o oca que el luido calien e de la placa de abajo asciende
po el cen o y el luido mas io desciende po las pa edes, p o ocando la dis ibución de empe a u as
mos ada en la igu a.
O o aspec o del p oblema de Rayleigh-Béna d que se puede es udia es modi ica las condiciones de
con o no en las pa edes.
Figu a( 6-7 )
Las pa edes la e ales amos a conside a las pe ec amen e conduc o as. La si uación es equi alen e a una
ca idad embebida en una plancha me álica con di e encia de empe a u a en e ca as. En la igu a 6-7 se
mues a el esquema gene al del p oblema.
El p oblema se esuel a idén icamen e al an e io a excepción de la condición de con o no en empe a u a
que en es e caso se educe a impone que la empe a u a asociada al mo imien o sea nula em las pa edes, de
es a o ma la empe a u a en las pa edes se á la de equilib io que se á la dis ibución que end á la plancha
me álica.
99
Figu a( 6-8 )
Obse amos como el e ec o de las pa edes conduc o as iene un e ec o amo iguado pa a la con ección ya
que acen úa más los e ec os de conducción de calo haciendo más di ícil la apa ición de un égimen
es aciona io. Cuando es e apa ece lo hace con es e pa ón di e en e al an e io debido a que la dis ibución de
empe a u as en las pa edes iene que se lineal po la condición de con o no de pa edes conduc o as que se
ha impues o.
Un es udio exhaus i o pa a de e mina el núme o de Rayleigh c í ico nos mos a ía que pa a una elación de
aspec o dada el núme o de Rayleigh c í ico pa a el caso de pa edes conduc o es se ía mayo que pa a pa edes
adiabá icas ya que como hemos comen ado es a condición de con o no acen úa los e ec os de conducción de
calo que amo iguan la ines abilidad p oducida po las ue zas de lo abilidad.
100
%==============ca idad con di e encia de empe a u as==============%
close all;clea all;
%=======================Pa áme os==========================%
nx=50;%in e alos de pa ición en x
ny=50;%in e alos de pa ición en y
n =150000 ;%in e alos de pa ición en
D=1; %ancho ca idad
D =0.0001 ; %inc emen o empo al
Dx=2*D/nx; %inc emen o eje x
Dy=D/ny; %inc emen o eje y
%=============NUMEROS ADIMENSIONALES==========================%
d =0.5; %INCREMENTO ADIMENSIONAL DE TEMPERATURA
M=0.1; %NÚMERO DE MACH
Re=30; %NÚMERO DE REYNOLDS
P =0.73; %NÚMERO DE PRANDTL
gamma=1.4; %COEFICIENTE DE EXPANSIÓN ADIABÁTICO
F =sq (d );%NÚMERO DE FROUDE
%============================TEMPERATURAS===============================
T2=1-d /2;
T1=1+d /2;
%============================COEFICIENTES===============================%
a1=D /Dx;
a2=D /Dy;
a3=D /F ^2;
a5=4*D /(3*Re*Dx^2);
a6=D /(Re*Dy^2);
a7=D /(Re*Dx^2);
a8=4*D /(3*Re*Dy^2);
a9=D /(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]=meshg id(x,y) ; %mallado
%======================CONDICIONES INICIALES===============%
P_u=10^-10*ones(nx+1,ny+1).*sin(2*pi*X).*cos(pi*Y);% elocidad esidual
di eccion x en los pun os del mallado
P_ =-10^-10*ones(nx+1,ny+1).*cos(2*pi*X).*sin(pi*Y);% elocidad esidual
di ección y en los pun os del mallado
P_Teq=(linspace(T1,T2,ny+1)'*ones(1,nx+1))'; %campo de empe a u a en cada
pun o del mallado equilib io
P_T=ze os(nx+1,ny+1); %campo de empe a u a DEL MOVIMIENTO en cada pun o
del mallado equilib io
P_ hoeq=P_Teq.^(gamma*M^2/d ^2-1);% densidad en los pun os del mallado
P_ ho=P_ hoeq;% densidad en los pun os del mallado
P_ hou=P_ ho.*P_u;% densidad*u en los pun os del mallado
P_ ho =P_ ho.*P_ ;%densidad* en los pun os del mallado
P_ hou2=P_ hou.*P_u;% densidad*u^2 en los pun os del mallado
P_ ho 2=P_ ho .*P_ ;%densidad* ^2 en los pun os del mallado
P_ hou =P_ hou.*P_ ;%densidad*u en los pun os del mallado
P_E=(P_T.*P_ ho-(P_ hoeq-P_ ho).*P_Teq)/gamma/(gamma-
1)/M^2+P_ hou2/2+P_ ho 2/2;
101
P_Eu=P_E.*P_u;
P_E =P_E.*P_ ;
u=10^-10*ones(nx+1,ny+1).*sin(2*pi*X).*cos(pi*Y);% elocidad esidual
di eccion x en los pun os del mallado
=-10^-10*ones(nx+1,ny+1).*cos(2*pi*X).*sin(pi*Y);% elocidad esidual
di ección y en los pun os del mallado
Teq=(linspace(T1,T2,ny+1)'*ones(1,nx+1))'; %campo de empe a u a en cada
pun o del mallado equilib io
T=ze os(nx+1,ny+1);
hoeq=Teq.^(gamma*M^2/d ^2-1);% densidad en los pun os del mallado
ho= hoeq;
hou= ho.*u;% densidad*u en los pun os del mallado
ho = ho.* ;%densidad* en los pun os del mallado
hou2= hou.*u;% densidad*u^2 en los pun os del mallado
ho 2= ho .* ;%densidad* ^2 en los pun os del mallado
hou = hou.* ;%densidad*u en los pun os del mallado
E=(T.* ho-( hoeq- ho).*Teq)/gamma/(gamma-1)/M^2+ hou2/2+ ho 2/2;
Eu=E.*u;
E =E.* ;
P_eq= hoeq.*Teq; %PRESIÓN EQUILIBRIO
Ra1=Re^2*P ; %NÚMERO RAYLEIGH STANDAR
Ra2=Re^2/F ^2*P *(max(max( ho))-min(min( ho))); %NÚMERO DE RAYLEIGH
CORREGIDO
pause
%========================BUCLE TEMPORAL===============================%
k=0;
con =1;
o q=1:n
k=k+1;
%===========================================PREDICTOR=======================
%
P_ ho(2:nx,2:ny)= ho(2:nx,2:ny)-a1*( hou(3:nx+1,2:ny)- hou(2:nx,2:ny))-
a2*( ho (2:nx,3:ny+1)- ho (2:nx,2:ny));
P_ hou(2:nx,2:ny)= hou(2:nx,2:ny)-a1*(gamma-1)*(E(3:nx+1,2:ny)-
E(2:nx,2:ny)- hou2(3:nx+1,2:ny)/2+ hou2(2:nx,2:ny)/2-
ho 2(3:nx+1,2:ny)/2+ ho 2(2:nx,2:ny)/2)-a1*( hou2(3:nx+1,2:ny)-
hou2(2:nx,2:ny))-a2*( hou (2:nx,3:ny+1)- hou (2:nx,2:ny))-
a10*u(2:nx,2:ny)+a5*(u(3:nx+1,2:ny)+u(1:nx-
1,2:ny))+a6*(u(2:nx,3:ny+1)+u(2:nx,1:ny-1))+a9*( (3:nx+1,3:ny+1)+ (1:nx-
1,1:ny-1)- (1:nx-1,3:ny+1)- (3:nx+1,1:ny-1));
P_ ho (2:nx,2:ny)= ho (2:nx,2:ny)-a2*(gamma-1)*(E(2:nx,3:ny+1)-
E(2:nx,2:ny)- hou2(2:nx,3:ny+1)/2+ hou2(2:nx,2:ny)/2-
ho 2(2:nx,3:ny+1)/2+ ho 2(2:nx,2:ny)/2)-a2*( ho 2(2:nx,3:ny+1)-
ho 2(2:nx,2:ny))-a1*( hou (3:nx+1,2:ny)- hou (2:nx,2:ny))-
a11* (2:nx,2:ny)+a7*( (3:nx+1,2:ny)+ (1:nx-
1,2:ny))+a8*( (2:nx,3:ny+1)+ (2:nx,1:ny-1))+a9*(u(3:nx+1,3:ny+1)+u(1:nx-
1,1:ny-1)-u(1:nx-1,3:ny+1)-u(3:nx+1,1:ny-1))-( ho(2:nx,2:ny)-
hoeq(2:nx,2:ny))*D /(F ^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*(gamma-
1)*0.5*(u(3:nx+1,2:ny).*( hou2(3:nx+1,2:ny)+ ho 2(3:nx+1,2:ny))-
u(2:nx,2:ny).*( hou2(2:nx,2:ny)+ ho 2(2:nx,2:ny)))+a1/(Dx*P *Re*M^2*(gamma-
1))*(T(3:nx+1,2:ny)-2*T(2:nx,2:ny)+T(1:nx-1,2:ny))-
a2*(gamma)*(E (2:nx,3:ny+1)-E (2:nx,2:ny))+a2*0.5*(gamma-

102
1)*( (2:nx,3:ny+1).*( hou2(2:nx,3:ny+1)+ ho 2(2:nx,3:ny+1))-
(2:nx,2:ny).*( hou2(2:nx,2:ny)+ ho 2(2:nx,2:ny)))+a2/(Dy*P *Re*M^2*(gamma-
1))*(T(2:nx,3:ny+1)-2*T(2:nx,2:ny)+T(2:nx,1:ny-1))-( ho(2:nx,2:ny)-
hoeq(2:nx,2:ny)*gamma/(gamma-1)).* (2:nx,2:ny)*D /(F ^2)-1/M^2/(gamma-
1)*P_eq(2:nx,2:ny).*(a1*(u(3:nx+1,2:ny)-u(2:nx,2:ny))+a2*( (2:nx,3:ny+1)-
(2:nx,2:ny)));
%===========================DECODIFICACIÓN=======================%
P_u=P_ hou./P_ ho;
P_ =P_ ho ./P_ ho;
Vmod2=P_ .^2+P_u.^2;
P_T=(P_ hoeq-P_ ho).*P_Teq./P_ ho+(P_E./P_ ho-Vmod2)*M^2*gamma*(gamma-1);
%===========================C.C===================================%
P_T(:,1)=0;%%Y=0 pa ed calien e
P_T(:,ny+1)=0;%%%%y=d pa ed ia
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_ ho(1,:)= ho(1,1:ny+1)-(a1/2)*(- hou(3,1:ny+1)+4* hou(2,1:ny+1)-
3* hou(1,1:ny+1));%%%%x=0 densidad
P_ ho(nx+1,:)= ho(nx+1,1:ny+1)+(a1/2)*(- hou(nx-
1,1:ny+1)+4* hou(nx,1:ny+1)-3* hou(nx+1,1:ny+1));%%%x=d densidad
P_ ho(2:nx,1)= ho(2:nx,1)-(a2/2)*(- ho (2:nx,3)+4* ho (2:nx,2)-
3* ho (2:nx,1));%%%%y=0 densidad
P_ ho(2:nx,ny+1)= ho(2:nx,ny+1)+a2/2*(- ho (2:nx,ny-1)+4* ho (2:nx,ny)-
3* ho (2:nx,ny+1));%%%y=d densidad
P_u(:,ny+1)=0;%%%% elocidad u y=d
P_ (:,ny+1)=0;%%%% elocidad y=d
P_u(:,1)=0;%%%% elocidad u y=0
P_ (:,1)=0;%%%% elocidad y=0
P_ (1,:)=0;
P_ (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_ hou=P_ ho.*P_u;% densidad*u en los pun os del mallado
P_ ho =P_ ho.*P_ ;%densidad* en los pun os del mallado
P_ hou2=P_ hou.*P_u;% densidad*u^2 en los pun os del mallado
P_ ho 2=P_ ho .*P_ ;%densidad* ^2 en los pun os del mallado
P_ hou =P_ hou.*P_ ;%densidad*u en los pun os del mallado
P_E=(P_T.*P_ ho-(P_ hoeq-P_ ho).*P_Teq)/gamma/(gamma-
1)/M^2+P_ hou2/2+P_ ho 2/2;
P_Eu=P_E.*P_u;
P_E =P_E.*P_ ;
%================================CORRECTOR================================%
103
ho(2:nx,2:ny)=0.5*(P_ ho(2:nx,2:ny)+ ho(2:nx,2:ny)-
a1*(P_ hou(2:nx,2:ny)-P_ hou(1:nx-1,2:ny))-a2*(P_ ho (2:nx,3:ny+1)-
P_ ho (2:nx,1:ny-1)));
hou(2:nx,2:ny)=0.5*( hou(2:nx,2:ny)+P_ hou(2:nx,2:ny)-a1*(gamma-
1)*(E(2:nx,2:ny)-E(1:nx-1,2:ny)- hou2(2:nx,2:ny)/2+ hou2(1:nx-1,2:ny)/2-
ho 2(2:nx,2:ny)/2+ ho 2(1:nx-1,2:ny)/2)-a1*(P_ hou2(2:nx,2:ny)-
P_ hou2(1:nx-1,2:ny))-a2*(P_ hou (2:nx,2:ny)-P_ hou (2:nx,1:ny-1))-
a10*P_u(2:nx,2:ny)+a5*(P_u(3:nx+1,2:ny)+P_u(1:nx-
1,2:ny))+a6*(P_u(2:nx,3:ny+1)+P_u(2:nx,1:ny-
1))+a9*(P_ (3:nx+1,3:ny+1)+P_ (1:nx-1,1:ny-1)-P_ (1:nx-1,3:ny+1)-
P_ (3:nx+1,1:ny-1)));
ho (2:nx,2:ny)=0.5*( ho (2:nx,2:ny)+P_ ho (2:nx,2:ny)-a2*(gamma-
1)*(E(2:nx,2:ny)-E(2:nx,1:ny-1)- hou2(2:nx,2:ny)/2+ hou2(2:nx,1:ny-1)/2-
ho 2(2:nx,2:ny)/2+ ho 2(2:nx,1:ny-1)/2)-a2*(P_ ho 2(2:nx,2:ny)-
P_ ho 2(2:nx,1:ny-1))-a1*(P_ hou (2:nx,2:ny)-P_ hou (1:nx-1,2:ny))-
a11*P_ (2:nx,2:ny)+a7*(P_ (3:nx+1,2:ny)+P_ (1:nx-
1,2:ny))+a8*(P_ (2:nx,3:ny+1)+P_ (2:nx,1:ny-
1))+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_ ho(2:nx,2:ny)-P_ hoeq(2:nx,2:ny))*D /(F ^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*(gamma-
1)*0.5*(P_u(2:nx,2:ny).*(P_ hou2(2:nx,2:ny)+P_ ho 2(2:nx,2:ny))-P_u(1:nx-
1,2:ny).*(P_ hou2(1:nx-1,2:ny)+P_ ho 2(1:nx-
1,2:ny)))+a1/(Dx*P *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_E (2:nx,2:ny)-
P_E (2:nx,1:ny-1))+a2*(gamma-
1)*0.5*(P_ (2:nx,2:ny).*(P_ hou2(2:nx,2:ny)+P_ ho 2(2:nx,2:ny))-
P_ (2:nx,1:ny-1).*(P_ hou2(2:nx,1:ny-1)+P_ ho 2(2:nx,1:ny-
1)))+a2/(Dy*P *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_ ho(2:nx,2:ny)-
P_ hoeq(2:nx,2:ny)*gamma/(gamma-1)).*P_ (2:nx,2:ny)*D /(F ^2)-1/M^2/(gamma-
1)*P_eq(2:nx,2:ny).*(a1*(P_u(2:nx,2:ny)-P_u(1:nx-
1,2:ny))+a2*(P_ (2:nx,2:ny)-P_ (2:nx,1:ny-1))));
%===========================DECODIFICACIÓN=======================%
u= hou./ ho;
= ho ./ ho;
Vmod2= .^2+u.^2;
T=( hoeq- ho).*Teq./ ho+(E./ ho-Vmod2)*M^2*gamma*(gamma-1);
%===========================C.C===================================%
T(:,1)=0;%%Y=0 pa ed calien e
T(:,ny+1)=0;%%%%y=d pa ed ia
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
ho(1,:)=(P_ ho(1,:)+ ho(1,:)-a1/2*(-P_ hou(3,:)+4*P_ hou(2,:)-
3*P_ hou(1,:)))*0.5;%%%%x=0
ho(nx+1,:)=(P_ ho(nx+1,:)+ ho(1,:)+a1/2*(-P_ hou(nx-
1,:)+4*P_ hou(nx,:)-3*P_ hou(nx+1,:)))*0.5;%%%x=d
ho(2:nx,1)=(P_ ho(2:nx,1)+ ho(2:nx,1)-a2/2*(-
P_ ho (2:nx,3)+4*P_ ho (2:nx,2)-3*P_ ho (2:nx,1)))*0.5;%%%%y=0
ho(2:nx,1+ny)=(P_ ho(2:nx,1+ny)+ ho(2:nx,1+ny)+a2/2*(-P_ ho (2:nx,ny-
1)+4*P_ ho (2:nx,ny)-3*P_ ho (2:nx,ny+1)))*0.5;%%%y=d
u(:,ny+1)=0;%%%% elocidad u y=d
(:,ny+1)=0;%%%% elocidad y=d
u(:,1)=0;%%%% elocidad u y=0
104
(:,1)=0;%%%% elocidad y=0
(1,:)=0;%%%%%X=0 d /dx=0
(nx+1,:)=0;%%%%%X=d d /dx=0
u(1,:)=0;%%%%%X=0 du/dx=0
u(nx+1,:)=0;%%%%%X=d du/dx=0
%====================================ACTUALIZACIÓN DE
VARIABLES===================%
hou= ho.*u;% densidad*u en los pun os del mallado
ho = ho.* ;%densidad* en los pun os del mallado
hou2= hou.*u;% densidad*u^2 en los pun os del mallado
ho 2= ho .* ;%densidad* ^2 en los pun os del mallado
hou = hou.* ;%densidad*u en los pun os del mallado
E=(T.* ho-( hoeq- ho).*Teq)/gamma/(gamma-1)/M^2+ hou2/2+ ho 2/2;
Eu=E.*u;
E =E.* ;
%==================REPRESENTACIÓN GRÁFICA Y CONVERGENCIA===================%
i k==500
D ;
q;
pel(:,:,con )=T;
con =con +1;
k=0;
end
H(q)=max(max(abs( )));
i max(max(u))>=20
q=n ;
pause
end
end
qui e (linspace(0,2,nx+1),linspace(0,1,ny+1),u'./sq ( '.^2+u'.^2), './sq (
'.^2+u'.^2))
105
7 CONCLUSIONES Y DESARROLLOS
FUTUROS
El obje i o p incipal de es e p oyec o ha sido demos a que, hoy en día, los es udian es de
disciplinas como la Dinámica de Fluidos y la T ans e encia de Calo pueden esol e
e icien emen e p oblemas ealis as haciendo uso de conocimien os básicos de p og amación y
de cálculo numé ico adqui idos en asigna u as cu sadas en años an e io es. El uso del
o denado pe mi e al alumno un mejo en endimien o de las leyes undamen ales de la
luidodinámica, sin que enga que ecu i a manipulaciones engo osas ni a simpli icaciones
que hagan pe de al p oblema odo su in e és, buscando una solución las ecuaciones que
gobie nan el p oblema.
Pa a ilus a es o se ha lle ado a cabo un análisis numé ico median e el mé odo de
MacCo mack aplicado a di e sos p oblemas luidodinámicos como el lujo al ededo de un
cilind o cuad ado que iaja a elocidad cons an e po un conduc o, la capa lími e é mica
al ededo del mismo y el égimen lamina no lineal de la con ección na u al de Rayleigh-
Béna d en el in e io de una ca idad ec angula .
En el capí ulo 3 se ha expues o el mé odo de MacCo mack explici o, el esquema mul i-e apa de
di e encias ini as, de segundo o den an o el iempo como en el espacio, que hemos u ilizado
pa a esol e los sis emas de ecuaciones que se han p esen ado en las di e en es aplicaciones.
Se han desa ollado dos ejemplos, el p ime o unidimensional, la esolución de la ecuación
lineal de onda, y el segundo, ya más en ocado a la esolución del abajo, un caso
bidimensional en el cual la con ección es a lle ada po la apade a, u iliza emos es os ejemplos
como alidado es del mé odo ya que como podemos obse a el acue do con los esul ados de
la bibliog a ía es excepcional.
En el capí ulo 4, se ha es udiado el lujo comp esible al ededo de un cilind o cuad ado,
p ime o ealizamos una compa ación de los coe icien es de sus en ación y esis encia con los
ob enidos en la li e a u a, ob eniendo un acue do excelen e. Después analizamos como a ec aba
la elación de aspec o del dominio compu acional a los esul ados pa a busca una solución de
comp omiso en e gas o compu acional y esul ados que se ap oximasen a un cilind o cuad ado
inme so en una co ien e inciden e. Una ez seleccionado el dominio compu acional
es udiamos el enómeno de ‘ o ex shedding’ y la ecuencia con la que los o bellinos se
emi ían en unción del núme o de Reynolds, el acue do con los da os de la bibliog a ía es muy
bueno eniendo en cuen a la limi ación compu acional que enemos. El inc emen o empo al
que nos exige la condición CFL es muy es ic i o y el espaciamien o de la malla se uel e
es ic i o cuando aumen amos el núme o de Reynolds lo que penaliza aún más la condición
CFL. Dos posibles soluciones pa a es o pod ían se un e inamien o de malla selec i o en zonas
donde sea necesa io, en es e caso en el en o no del cuad ado se pod ía hace una malla más
ina, y en el en o no de las secciones de en ada y salida y pa edes la e ales se pod ía ene una
malla más g uesa ya que lo único que que emos en esos pun os es impone condiciones de
con o no y el desa ollo del mé odo de MacCo mack modi icado que se explica b e emen e en
el capí ulo es que elaja la condición CFL al di idi el p oceso en e apas unidimensionales.
En el capí ulo 5, se ha desa ollado una b e e aplicación como p ime con ac o a la inclusión de
la empe a u a en las ecuaciones. Los pe iles de empe a u as ob enidos e an consecuen es con
lo espe ado y la es imación de ó denes de magni ud del espeso de la capa lími e é mica
ambién lo e an. Es e p oblema equie e una capacidad compu acional ele ada po eso no se