TRABAJO DE FIN DE GRADO - Archivo Digital UPM -...

Escuela Técnica Superior de Ingeniería Industrial

Universidad Politécnica de Madrid

Grado en Ingeniería en Tecnologías Industriales

TRABAJO DE FIN DE GRADO

TEORÍA DEL FLUJO POTENCIAL Y CAPA

LÍMITE PARA LA CARACTERIZACIÓN DEL

MOVIMIENTO FLUIDO CON APLICACIÓN

EN ESTENOSIS ARTERIAL

Marta Herrero Hernanz

Junio 2016

Agradecimientos

A mis padres y a mis hermanos, por brindarme todo su interés y conanza.

A Leo, por estar ahí en los momentos más necesarios.

A mis amigas, futuras médicos, por su disposición y ayuda otorgadas.

A Javier y a Jaime, por su incondicional apoyo, dedicación y ayuda.

Resumen

El presente proyecto tiene como objetivo el estudio y caracterización del movimientouido en determinadas situaciones en las que los efectos viscosos del uido sean despreciablesen la mayor parte del dominio y estén reducidos a una región estrecha del mismo. Además elnúmero de Reynolds caracteriza el régimen del movimiento, pudiendo ser laminar o turbulentoen función del valor que tome. Este número adimensional será de vital importancia a lo largodel proyecto, además de otros parámetros uidodinámicos como el esfuerzo cortante o elcoeciente de fricción.

Las principales ecuaciones que rigen la Mecánica de Fluidos son las ecuaciones de con-servación de masa, de cantidad de movimiento y de energía. Su resolución directa es compleja,por lo que frecuentemente para acelerar el proceso del cálculo se realizan simplicaciones delmodelo. Para el caso de este proyecto, se abordarán ujos con densidad constante, estacio-narios y con fuerzas másicas despreciables. La primera simplicación que aparece es que losproblemas mecánico y térmico están desacoplados, por lo que solo es necesario resolver elproblema mecánico donde el campo de velocidades y de presiones está dado por las famosasecuaciones de Navier Stokes:

∇·~v = 0

ρ(~v·∇)~v = −∇p+ µ4~v

El tipo de movimiento que se pretende estudiar puede observarse grácamente en lagura 1. Existe una primera región cerca del cuerpo donde los esfuerzos viscosos son predo-minantes, conocida como capa límite, mientras que fuera de ésta los esfuerzos viscosos sondespreciables. Esta distinción es la división natural del proyecto, pues en la zona de la capalímite se aplicará la teoría de la capa límite laminar mientras que en la zona exterior a éstase aplicará teoría del ujo potencial, un caso particular de movimientos ideales. Además elestudio de la capa límite es de gran interés, pues es la responsable de la resistencia a lafricción que ofrecen todos los cuerpos inmersos en un uido. El desprendimiento de la capalímite es un fenómeno asociado a la misma que ocurre si ésta está sometida a gradientes depresión adversos, y que origina vórtices que conllevan disipación de energía.

La teoría del ujo potencial, el despreciar los efectos viscosos, requerirá que los númerosde Reynolds sean altos. Por otro lado, la teoría de la capa límite laminar ha de aplicarse aujos con números de Reynolds que se incluyan en el régimen laminar. Por lo que el principalobjetivo del proyecto es crear un modelo matemático que permita el acoplamiento de ambas

Figura 1: Flujo alrededor de un cuerpo

regiones para poder así conocer el campo uido completo alrededor de un cuerpo. El númerode Reynolds del uido ha de ser lo sucientemente alto como para poder despreciar los efectosviscosos, si bien que no sobrepase el rango de régimen laminar.

Las ecuaciones de Navier Stokes se adaptan a cada una de las zonas, por lo que esnecesario el tratamiento matemático previo de las mismas. La simplicación de las ecuacionessegún la teoría del ujo potencial permitirá obtener unas ecuaciones en términos de la funciónde corriente y de la función potencial, y a partir de éstas se obtienen los campos de presiones yvelocidades en la zona exterior a la capa límite. En cambio, la teoría de la capa límite laminarestá regida por unas ecuaciones cuyas incógnitas son las velocidades del uido horizontal yvertical. El acoplamiento entre ambas soluciones es posible gracias a que la presión se conservaen dirección perpendicular a la supercie del cuerpo.

Para poder llevar a cabo la resolución y acoplamiento de las teorías se procede a laresolución numérica de las ecuaciones ya que la resolución analítica no es posible. Además seadimensionalizan todas las variables para mayor comodidad en su tratamiento. Primeramentese resolverá la parte de ujo potencial, para lo cual se ha utilizado el programa FreeFem++,que resuelve las ecuaciones de la teoría del ujo potencial mediante elementos nitos. Eneste caso el mallado lo crea el propio programa, siendo solo necesaria la denición de lageometría y de las condiciones de contorno sobre las funciones de corriente. La resolución dacomo resultado los campos de presiones y velocidades, y como se ha dicho, será a través dela presión el acoplamiento de las teorías.

En cuanto a la teoría de la capa límite laminar, las ecuaciones requieren un mayortratamiento previo a su resolución. Será necesario discretizar las ecuaciones adimensiona-lizadas, lo cual se realizará por diferencias nitas, además de denir el mallado. Medianteel renamiento del mismo se controla la precisión del método y, consecuentemente, el errorcometido. Las condiciones de contorno en la zona de la capa límite serán la condición de nodeslizamiento en la supercie del cuerpo y la velocidad y presiones en la zona exterior a lacapa límite. La resolución de estas ecuaciones se realiza en Matlab.

Una vez planteado el modelo matemático y programados los métodos numéricos parasu resolución se procede a su vericación. Para ello se resuelve el ujo alrededor de una placaplana y alrededor de un cilindro. Se han elegido estas geometrías porque se dispone de lasolución analítica de las ecuaciones de la capa límite laminar con condiciones de contorno deujo potencial en la bibliografía consultada [12].

Para el caso de la placa plana, la teoría potencial es trivial, pues la presión es la mismaa lo largo de todo el campo uido. Tras la resolución de la capa límite laminar se obtienenlos perles de velocidades en todos los puntos de la placa plana. Para vericar la correctaresolución de este caso, como se ha dicho, se comparará con la solución pseudoanalítica.La comparación de los perles de velocidades y del coeciente de fricción medio permitenconcluir que se ha resuelto con exactitud y exitosamente este caso. Es de esperar que no hayadesprendimiento de la capa límite al no existir gradientes de presión adversos.

Análogamente se resuelve el caso de ujo alrededor de un cilindro, con la misma di-námica, variando la condiciones de contorno, adaptándolas a la geometría de este problema.En este caso la solución de la teoría del ujo potencial no es trivial, pues las presiones va-rían a lo largo de la supercie del cilindro y existen gradientes adversos, por lo que se prevéque la capa límite se desprenda. Tras la obtención del campo de velocidades es posible lacomparación con los perles de velocidades teóricos así como el punto de desprendimiento.Éste último tiene un error del 3,05 % con respecto al esperado, obteniéndose un ángulo dedesprendimiento de 105 medido desde el punto de remanso.

Además se ha probado un método de control de la capa límite: la succión homogéneaa lo largo del cuerpo. Este método modica la capa límite de forma articial, y consiste enla absorción de partículas uidas a través de una de las paredes del cuerpo sobre el que seforma la capa límite. De esta manera se reduce el espesor de la capa límite, por lo que seretrasa su crecimiento y posible desprendimiento como se muestra en la gura 2.

Figura 2: Succión de la capa límite [12]

Se ha probado tanto en la placa plana como en el cilindro. Para el primero de ellos sedispone de solución analítica con la que comparar, mientras que en el cilindro se ha realizadouna comparación cualitativa con el caso de succión nula. A modo de ejemplo se presentan losresultados obtenidos para el caso de una placa plana en la gura 3, donde pueden verse losperles de velocidades en distintos puntos para una placa de longitud unidad. Como puedeapreciarse la velocidad es nula en la supercie de la placa, y la capa límite crece conforme seavanza en longitud en la placa. Comparando los perles obtenidos con y sin succión puedeapreciarse un claro retraso en el crecimiento de la capa límite debido a la succión.

Tras resolver las geometrías de placa plana y cilindro y vericar la validez del método,se ha procedido al desarrollo de la aplicación nal del mismo. En los últimos años, la impor-tancia de la ingeniería en la medicina ha ido creciendo hasta ser, a día de hoy, dos ciencias

0 0.5 1

u

0

1

2

3

4

5

6

7

Altu

ra

x=0.25

0 0.5 1

u

x=0.5

0 0.5 1

u

x=0.75

0 0.5 1

u

x=1

Sin succiónCon succión

Figura 3: Perles de velocidades en distintos puntos de la placa plana

estrechamente relacionadas. Es por ello que se ha elegido una aplicación médica: la estenosisarterial. Se trata de una enfermedad que afecta a los vasos sanguíneos de tal manera que cier-tas arterias se estrechan como se aprecia en la gura 4. El objetivo de esta parte es conocerel ujo alrededor de la estenosis, para lo cual el método desarrollado ha de modicarse enciertos aspectos. Esto se debe a que en el caso del ujo sanguíneo los efectos viscosos no sondespreciables, por lo que la teoría del ujo potencial no es la más adecuada para imponer enla zona exterior a la capa límite. Por otro lado, el estrechamiento provoca la aceleración delujo y localmente los efectos de inercia predominan sobre los viscosos, por ello cabe hablarde la capa límite en el estrechamiento del vaso y no antes de la estenosis, donde el perl develocidades el de tipo Poiseuille como se puede ver en la gura 4.

Figura 4: Geometría real y aproximada de la estenosis arterial

Para vericar la validez de lo que se realiza en esta parte del proyecto, se dispone dedatos experimentales sobre el esfuerzo cortante máximo a lo largo de la estenosis, con los quese comparará. Se resuelven cinco casos distintos, siendo todos ellos variaciones semejantessobre una misma geometría. Además se realiza una corrección sobre el perl de velocidadesen la garganta de la estenosis, para lo cual se explica el concepto de espesor de desplazamiento.La tensión máxima sobre la estenosis se muestra para tres de los casos en la gura 5.

102 103

Re

100

101

102

103

τ e

n di

nas/

cm2

τ máximo en función del Reynolds

Upstream of throat M3



Experimental

Figura 5: Resultados obtenidos: esfuerzo cortante máximo en la estenosis

Por último, para relacionar el método utilizado en este apartado con el desarrolladoen el proyecto se resuelven estos casos para el número de Reynolds que permite seguir en elrégimen laminar para el caso del ujo sanguíneo. Para este caso la aplicación de la teoría deujo potencial sí es correcta, y la comparación de los resultados obtenidos con ésta y conla teoría de Poiseuille se realiza a través de los puntos de desprendimiento. Éstos estaránmás adelantados para la resolución con condiciones de contorno de ujo potencial, pues sedesprecian los efectos viscosos en el exterior de la capa límite y por ello la adherencia de lamisma a la pared del vaso es menor, y se desprende antes.

Finalmente, se presentan una serie de conclusiones generales que justican la validezdel método desarrollado. Una de las principales ventajas del método desarrollado es el escasotiempo de cálculo de los programas de ordenador junto con la exactitud del método, aceptablepara las aplicaciones que se han tratado. La conciliación entre las teorías de capa límitelaminar y ujo potencial dan lugar a una herramienta de cálculo de gran aplicación, si bienes necesario conocer sus limitaciones.

Palabras clave: Reynolds, función de corriente, capa límite laminar, ujo potencial,esfuerzo cortante estenosis.

Códigos UNESCO220403 - ujo de uidos220504 - mecánica de uidos250121 - simulación numérica120600 - análisis numérico330809 - ingeniería sanitaria320702 - arteriosclerosis

Índice general

1. Introducción 1

1.1. Antecedentes y motivación . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2. Conceptos básicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2.1. Descripción cuantitativa del movimiento de un uido alrededor de uncuerpo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2.2. Descripción cualitativa del movimiento de un uido alrededor de uncuerpo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.3. Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.4. Metodología . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.5. Estructura del proyecto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2. Teoría del ujo potencial 9

2.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.2. Condiciones de contorno . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.3. Paradoja de D'Alambert . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3. Teoría de la capa límite laminar 15

3.1. Simplicación de las ecuaciones por análisis dimensional . . . . . . . . . . . 15

3.2. Condiciones de contorno . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3.3. Desprendimiento de la capa límite . . . . . . . . . . . . . . . . . . . . . . . . 17

3.4. Espesor de la capa límite para la placa plana: ecuación de Karman . . . . . . 18

i

ÍNDICE GENERAL

3.5. Coeciente de fricción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.6. Solución analítica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.6.1. Caso: placa plana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.6.2. Solución pseudoanalítica . . . . . . . . . . . . . . . . . . . . . . . . . 22

4. Modelo matemático 27

4.1. Introducción a los métodos numéricos . . . . . . . . . . . . . . . . . . . . . . 27

4.2. Adimensionalización de las variables . . . . . . . . . . . . . . . . . . . . . . . 27

4.3. Datos del problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

4.4. Campo de velocidades según la teoría del ujo potencial. . . . . . . . . . . . 29

4.4.1. Mallado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

4.4.2. Discretización y datos . . . . . . . . . . . . . . . . . . . . . . . . . . 30

4.4.3. Resolución . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

4.5. Campo de velocidades en la capa límite laminar . . . . . . . . . . . . . . . . 31

4.5.1. Mallado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

4.5.2. Discretización y datos . . . . . . . . . . . . . . . . . . . . . . . . . . 33

4.5.3. Resolución . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

4.5.4. Coeciente de fricción . . . . . . . . . . . . . . . . . . . . . . . . . . 38

5. Flujo alrededor de una placa plana 39

5.1. Flujo potencial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

5.2. Capa límite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

5.2.1. Comparación con la solución pseudoanalítica . . . . . . . . . . . . . . 43

5.2.2. Coeciente de fricción . . . . . . . . . . . . . . . . . . . . . . . . . . 44

5.3. Succión homogénea en la placa plana . . . . . . . . . . . . . . . . . . . . . . 47

5.3.1. Solución analítica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

5.3.2. Flujo potencial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

5.3.3. Capa límite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

ii ETSII-UPM

TEORÍA DEL FLUJO POTENCIAL Y CAPA LÍMITE

6. Flujo alrededor de un cilindro 55

6.1. Características del ujo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

6.2. Flujo potencial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

6.3. Capa límite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

6.4. Succión homogénea en el cilindro . . . . . . . . . . . . . . . . . . . . . . . . 64

6.4.1. Flujo potencial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

6.4.2. Capa límite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

7. Aplicación: ujo en arterias con estenosis 69

7.1. Geometría . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

7.2. Características del ujo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

7.3. Datos experimentales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

7.4. Metodología . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

7.5. Resolución del ujo sanguíneo con condiciones de contorno de ujo de Poiseuille 75

7.5.1. Corrección de la velocidad con el espesor de desplazamiento . . . . . 75

7.5.2. Datos del problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

7.5.3. Cálculo del esfuerzo cortante . . . . . . . . . . . . . . . . . . . . . . . 80

7.5.4. Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

7.5.5. Comprobación de hipótesis . . . . . . . . . . . . . . . . . . . . . . . . 81

7.6. Extrapolación a mayores números de Reynolds: introducción de la teoría delujo potencial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

7.6.1. Flujo potencial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

7.6.2. Capa límite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

7.7. Vericación del método . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

7.8. Conclusiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

7.9. Mejoras en la geometría . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

Marta Herrero Hernanz iii

ÍNDICE GENERAL

Conclusiones 90

A. Planicación temporal y presupuesto 93

A.1. Planicación temporal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

A.1.1. Estructura de Descomposición del Proyecto . . . . . . . . . . . . . . 93

A.1.2. Diagrama de Gantt . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

A.2. Presupuesto del proyecto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

B. Reexión general sobre aspectos de responsabilidad social del trabajo 99

C. Código fuente 101

C.1. Placa plana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

C.1.1. Solución pseudoanalítica . . . . . . . . . . . . . . . . . . . . . . . . . 101

C.1.2. Solución por métodos numéricos . . . . . . . . . . . . . . . . . . . . . 102

C.2. Cilindro . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

C.2.1. Sin succión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

C.2.2. Con succión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

C.3. Estenosis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

C.3.1. Capa límite con condiciones de contorno de Poiseuille . . . . . . . . . 113

C.3.2. Capa límite con condiciones de contorno de ujo potencial . . . . . . 117

Índice de guras 121

Índice de cuadros 124

Nomenclatura 126

Bibliografía 129

iv ETSII-UPM

Capítulo 1

Introducción

1.1. Antecedentes y motivación

Las aplicaciones de la Mecánica de Fluidos a la vida real son de gran importancia yabarca ámbitos tan distintos como el diseño de vehículos o la medicina. Fenómenos como eldesprendimiento de la capa límite, la sustentación de una aeronave o las ondas de choquehacen interesante el estudio de los uidos.

El ambicioso objetivo de conocer el comportamiento de un uido alrededor de un cuer-po, ya sea éste un vehículo, los álabes de un turbocompresor, o cualquier otro, es alcanzablesi se resuelven las ecuaciones de Navier Stokes. Estas ecuaciones describen el movimiento deun uido en un entorno, si bien su resolución no es sencilla, por lo que en gran parte de lasaplicaciones se requieren simplicaciones matemáticas previas a su tratamiento.

La fenomenología que rodea a la capa límite es de gran interés, pues es la responsablede la resistencia a la fricción que ofrecen todos los cuerpos inmersos en un uido. La capalímite fue planteada en 1904 por Prandtl, al que se le conoce como el Padre de la Mecánica deFluidos Moderna. La teoría de la capa límite supone la conciliación de la hidrodinámica y lahidráulica, pues permite entender cómo una pequeña región viscosa afecta a las característicasglobales del ujo.

Por otro lado, la aplicación de las leyes que rigen los uidos al caso de la sangre permiteun acercamiento más o menos realista a la fenomenología sanguínea. En los últimos años laaplicación de la Mecánica de Fluidos en la medicina ha ido creciendo, lo cual es un granapoyo para la hemodinámica y todo ello se logra a través de la simulación.

1.2. Conceptos básicos

Un uido es una sustancia que se deforma de forma continua cuando está sometidoa un esfuerzo cortante. De cara a un estudio de ingeniería, cabe comentar la importancia

1

CAPÍTULO 1. INTRODUCCIÓN

del régimen de movimiento de un uido, que puede ser laminar o turbulento. El primerode ellos alude a un movimiento ordenado, con una variación progresiva y gradual de laspropiedades uidas, mientras que el segundo de ellos describe el movimiento de un uido deforma desordenada, formando torbellinos de gran variedad de tamaños. La caracterización delrégimen de movimiento viene dada por el número de Reynolds (1883), que es adimensionaly se dene como:

Re =ρvD

µ

donde ρ es la densidad del uido, v la velocidad del uido, D la longitud característicadel problema y µ la viscosidad dinámica del uido.

De esta manera, a modo de ejemplo, para el caso de una tubería circular para Re < 2000se tiene ujo laminar, y cualquier turbulencia que se produzca es eliminada por la viscosidad,mientras que para Re > 3000 el régimen es turbulento. Entre estos dos valores se da lazona de transición entre ambos números. El límite entre laminar y turbulento varía de unageometría a otra.

1.2.1. Descripción cuantitativa del movimiento de un uido alrede-

dor de un cuerpo

A continuación se presentan las ecuaciones que permiten describir y cuanticar el movi-miento de los uidos, a través de la conservación de masa, cantidad de movimiento y energía,que constituyen un sistema de cinco ecuaciones con cinco incógnitas.

Conservación de la masa

Dρ

Dt+ ρ∇·~v = 0

Conservación de la cantidad de movimiento

ρD~v

Dt= ∇·τ + ρ~fm

donde ~fm son las fuerzas másicas y τ el tensor de esfuerzos viscosos.

Conservación de la energía

ρDe

Dt= −p∇·~v + Φv −∇~q +Qr +Qq

2 ETSII-UPM


donde e es la energía interna, Φv la función de disipación de Rayleigh, ~q el calor porconducción y Qr y Qq los calores de reacción y químicos.

Frecuentemente la densidad del uido es constante a lo largo de la aplicación, es decir,el uido es incompresible. Imponiendo esta condición se obtienen las conocidas ecuaciones deNavier Stokes:

∇·~v = 0

ρD~v

Dt= −∇p+ µ4~v + ρ~fm

ρDe

Dt= Φv −∇~q +Qr +Qq

En el caso de líquidos, donde ρ y µ se pueden considerar constantes y uniformes enel movimiento uido, los problemas mecánico y térmico están desacoplados. Es por ello quees posible resolver el problema mecánico dejando de lado el problema térmico, por lo que setrabajará con las dos primeras ecuaciones. Además la existencia o inexistencia de evolucióntemporal permitirá diferenciar entre ujos transitorios o estacionarios. En este proyecto seconsiderarán todos los procesos estacionarios, es decir, el valor de las variables uidas sonconstantes a lo largo del tiempo. Además de suponer las fuerzas másicas despreciables. Así,las ecuaciones con las que se va a trabajar son las siguientes:

∇·~v = 0 (1.1)

ρ(~v·∇)~v = −∇p+ µ4~v (1.2)

En el caso de movimiento uido alrededor de cuerpos, un resultado interesante de medires el coeciente de arrastre y sustentación. De manera general, la resistencia y la sustentaciónson las resultantes de las fuerzas de presión Fp y las debidas a esfuerzos cortantes (τ) o fuerzasde fricción Ff sobre la supercie del cuerpo. Por lo tanto, si se determinan las distribucionesde presiones y de esfuerzos cortantes sobre la supercie, se pueden obtener por integraciónlas fuerzas de presión y de fricción.

−→Ff =

ˆA

(p− p∞)−→n dA

−→Fp =

ˆA

τ−→t dA

donde −→n y−→t son los vectores normal y tangencial respectivamente a la supercie.

Cada una de estas fuerzas se puede dividir en dos, una en sentido normal al ujoincidente y otra en sentido paralelo al mismo, lo cual se reeja en la gura 1.1 para el casode un perl aerodinámico. Así, la fuerza en sentido vertical se llama fuerza de sustentación(L) o lift , y la resultante en el sentido horizontal es la fuerza de arrastre (D) o drag.

Marta Herrero Hernanz 3


Figura 1.1: Fuerzas de fricción y de presión sobre un perl aerodinámico

Por comodidad, se denen los coecientes de fricción Cf y de presión Cp, que son estasfuerzas adimensionalizadas:

Cf =Ff

12ρU2∞A

Cp =Fp

12ρU2∞A

De igual manera se pueden denir los coecientes de arrastre CD y de sustentación CLcomo:

CD =D

12ρU2∞A

CL =L

12ρU2∞A

donde A es el área proyectada sobre el plano horizontal, y Ff y Fp las fuerzas de presióny de fricción en sentido horizontal.

1.2.2. Descripción cualitativa del movimiento de un uido alrededor

de un cuerpo

Las ecuaciones de Navier Stokes simplicadas permiten conocer el campo de velocidadesde un uido y presiones como se ha dicho, si bien su resolución no es sencilla, ya que sonecuaciones no lineales en derivadas parciales. En pocos casos, realizando simplicaciones enlas mismas, se pueden obtener los campos de velocidad y presiones de forma analítica.

4 ETSII-UPM


Sea un cuerpo sumergido en un uido con velocidad, las partículas de éste rodean elcuerpo y siguen su trayectoria de tal manera que recuperan su velocidad inicial tras alejarsedel cuerpo. Sin embargo, alrededor del cuerpo el campo de velocidades es alterado, y elestudio de este fenómeno es de especial interés.

Forzosamente, la velocidad del uido debe ser nula en el contacto con el cuerpo, y nonula en una región más o menos próxima al contorno del cuerpo, como puede verse en lagura 1.2. La región de transición se conoce como capa límite, en la que los esfuerzos viscososson predominantes. Es importante tener en cuenta que la capa límite solo existe si el númerode Reynolds es sucientemente alto. En resumen, este proyecto tiene como objetivo el estudiodel campo de velocidades y presiones en capa límite y fuera de ella, mediante la simplicaciónde las ecuaciones de Navier Stokes.

Así, en una región alejada del cuerpo los esfuerzos viscosos son despreciables, y elanálisis se realizará a través de la teoría de ujo potencial. En cambio, en la capa límitehabrá que tenerlos en cuenta.

Figura 1.2: Flujo alrededor de un cuerpo

1.3. Objetivos

Este Trabajo de Fin de Grado pretende estudiar las simplicación matemática de lasecuaciones de Navier Stokes para su resolución, requiriendo conocimientos de Mecánica deFluidos, Matemáticas y programación de código fuente. Los objetivos claros de este proyectoson los siguientes:

Simplicación de las ecuaciones de Navier Stokes mediante la división del estudio delujo en dos partes.

Programación de las ecuaciones simplicadas en código fuente para obtener los camposde velocidades y de presiones del uido.

Resolución de las ecuaciones simplicadas y análisis de los resultados obtenidos.

Vericación del método desarrollado a través de la resolución de dos casos de geometríaconocida.



Estudio de fenomenología asociada a la capa límite laminar. Control de la capa límite:succión.

Aplicación del modelo matemático al caso de ujo sanguíneo para el estudio de laestenosis arterial. Comprobación de hipótesis.

Aprendizaje en profundidad de los conceptos con los que se va a trabajar mediante elestudio de los uidos y sus ecuaciones asociadas, así como gracias al uso de métodosnuméricos para la simplicación de las ecuaciones.

Con todo ello será posible, nalmente, un mejor entendimiento del comportamiento de losuidos en determinados escenarios y la importancia de su estudio. El objetivo más generales, por lo tanto, vericar la validez del modelo simplicado matemático a través de casosconocidos de la literatura y aplicarlo al caso de ujo en vasos sanguíneos.

1.4. Metodología

Tras haber expuesto los objetivos claros del proyecto, se procede a la descripción decómo se van a lograr.

Como se ha comentado, el estudio del ujo se divide en dos partes diferenciadas enbase a las características del uido en cada una de ellas. Alrededor del cuerpo los efectosviscosos son predominantes, por lo que las ecuaciones de Navier Stokes serán simplicadasteniendo esto en cuenta. Esta zona es la región donde existe la capa límite, y se aplicaráteoría de la capa límite laminar. El hecho de que sea laminar requiere que los números deReynolds se correspondan a la región de ujo laminar, pues si fueran de ujo turbulentolas ecuaciones no serían válidas. Por otro lado, en la zona exterior de la misma y, por lotanto, más alejada del cuerpo, los efectos viscosos son despreciables, por lo se aplicación dela teoría del ujo potencial está justicada como se explica en el capítulo correspondiente.En esta región los números de Reynolds son altos, tendiendo a innito, debido al despreciode los efectos viscosos. Por lo tanto, será necesaria la elección de un número de Reynoldssucientemente alto para poder aplicar teoría del ujo potencial si bien que se encuentredentro del régimen laminar.

Una vez simplicadas las ecuaciones matemáticas se procede a su resolución, para locual se han utilizado los programas Matlab y FreeFem++. En primer lugar, las ecuacioneshan de ser adaptadas a los códigos fuente de estos programas, para lo cual es necesariauna primera discretización de las mismas y, para mayor comodidad, su adimensionalización.Además ha de denirse el mallado para ambos programas.

FreeFem++ se utiliza para resolver las ecuaciones de la teoría del ujo potencial, y sebasa en resolución por elementos nitos. Una vez obtenidos los resultados se introducen éstoscomo condición de contorno de las ecuaciones de la capa límite laminar. Se resolverán las

6 ETSII-UPM


ecuaciones mediante diferencias nitas utilizando Matlab, por lo que son necesarios conoci-mientos de métodos numéricos para su correcto planteamiento. El acoplamiento entre ambaszonas se realiza, como se ha dicho, a través de la imposición de las condiciones de la capalímite laminar, lo cual es posible gracias a la conservación de la presión en dirección normala la supercie del cuerpo.

Conocido el modelo matemático, se aplica éste a los casos de ujo alrededor de unaplaca plana y ujo alrededor de un cilindro. Se eligen estas dos geometrías dado que se disponede la solución analítica de estas ecuaciones, por lo que será posible la comparación de losresultados obtenidos. Esta comparación se realizará a través de los perles de velocidades yel cálculo de curvas que aparecen en la literatura consultada [12].

Además se presenta un método de control de la capa límite: la succión. El estudio deeste caso requiere de la adaptación de determinadas condiciones de contorno como se explicaen el contenido del proyecto. Igualmente, se utiliza el modelo matemático desarrollado parala obtención del campo de velocidades. Las diferencias entre los casos de succión homogéneaa lo largo del cuerpo y sin succión se llevan a cabo a través de la comparación de los perlesde velocidades y de los puntos de desprendimiento de la capa límite.

Una vez vericada la validez del modelo matemático gracias a la resolución de casossencillos se da el salto a una aplicación más realista y de gran interés: ujo en arterias conestenosis. Para el estudio del uido alrededor de este tipo de geometrías se han realizadoalgunos cambios en el modelo matemático. En primer lugar, al tratarse de ujo sanguíneo,los números de Reynolds son bajos, por lo que la aplicación de la teoría potencial en la zonaexterior a la capa límite no es la mejor opción. El estrechamiento provoca la aceleracióndel ujo y localmente los efectos de inercia predominan sobre los viscosos, por ello cabehablar de la capa límite en el estrechamiento del vaso. En cambio, en la parte del vaso sinestrechamiento el ujo sanguíneo se considera del tipo de Poiseuille, por lo que en esta zonano se puede hablar de capa límite. En este apartado, consecuentemente, solo será necesarioel código fuente desarrollado en Matlab al no ser necesaria la resolución del campo de ujopotencial.

Se dispone de datos experimentales de cinco casos de avance de la enfermedad paranúmeros de Reynolds de 100, 500 y 1000. Por lo tanto, la resolución del ujo sanguíneo enarterias con estenosis se ha realizado para estas cinco geometrías, siendo todas ellas seme-jantes. Además se han realizado hipótesis sobre el carácter del perl de velocidades en lagarganta de la estenosis, las cuales se comprueban tras su resolución.

Finalmente, se ha realizado una comparación entre los resultados que se obtienen impo-niendo condiciones de contorno de ujo de Poiseuille y los que resultan de imponer condicionesde contorno de ujo potencial. Este comparación se lleva a cabo en los cinco estudiados, sibien el número de Reynolds se aumenta hasta 2000 para que sea válida la aplicación dela teoría potencial, siendo este el valor límite del régimen laminar en vasos sanguíneos. Losmétodos se comparan a través del estudio de los puntos de desprendimiento de la capa límite.



1.5. Estructura del proyecto

Tras conocer la temática y objetivos del proyecto, se procede a la presentación de laestructura del mismo, la cual permitirá claricar cómo se llevará a cabo el estudio.

En el Capítulo I se ha realizado una breve introducción a los uidos, y se ha expli-cado tanto cualitativa como cuantitativamente cómo es el ujo alrededor de un cuerpo. Sepresentan las ecuaciones de Navier Stokes con las que se va a trabajar a lo largo de todo eldesarrollo, así como conceptos básicos relacionados con los uidos.

En el Capítulo II se desarrolla la teoría del ujo potencial en profundidad, explicándosela procedencia de las ecuaciones que la gobiernan.

La teoría de la capa límite laminar es desarrollada en el Capítulo III, explicándose lasimplicación de las ecuaciones por análisis dimensional y el modo de acoplamiento con lateoría del ujo potencial. También incluye la explicación de parámetros relacionados con lacapa límite laminar de gran interés como el coeciente de fricción o el esfuerzo cortante. Porúltimo se explica la resolución analítica de las ecuaciones de la capa límite para el caso de laplaca plana, que se utilizará posteriormente para comparar con los resultados obtenidos conel modelo matemático desarrollado.

Es en el Capítulo IV donde se explica dicho modelo matemático. Se realiza una breveintroducción a los métodos numéricos y a continuación se procede a la adimensionalizaciónde las variables y la presentación de los datos del problema. En este capítulo se explica cómoFreeFem++ y Matlab resuelven las ecuaciones, para lo cual se habla del mallado de ambosproblemas y de la discretización y resolución de las ecuaciones.

Una vez presentado el modelo matemático se da pie a la resolución de los casos degeometría sencilla comentados. En el Capítulo V se desarrolla el caso de ujo alrededor deuna placa plana, explicándose las partes de ujo potencial y capa límite y comparando losresultados obtenidos con los de la solución analítica. Por último se presenta la succión enesta geometría. El Capítulo VI tiene una estructura similar al anterior si bien en este caso seresuelve el ujo alrededor de un cilindro.

Por último, en el Capítulo VII se desarrolla toda la parte de aplicación del modelomatemático. Se realiza una introducción a la estenosis arterial y se presenta la geometríay datos del problema. A continuación se explican las diferencias entre el método de ujoexterior a la capa límite de tipo Poiseuille y ujo exterior de tipo potencial, así como lametodología para abordar su resolución. Tras ésta se explica la extrapolación del método anúmeros de Reynolds altos.

Por último, se presentan las conclusiones generales del proyecto y las posibles lineasfuturas del mismo. En los Anexos se incluye la planicación temporal y presupuesto delproyecto, además de una breve reexión sobre el impacto social del trabajo. También seincluyen los principales códigos fuente utilizados para la resolución de los problemas.

8 ETSII-UPM

Capítulo 2

Teoría del ujo potencial

2.1. Introducción

Las ecuaciones de Navier Stokes permiten conocer el campo de velocidades y presionesde un uido en una región del espacio. Su resolución es complicada, por lo que habrá quehacer simplicaciones para ello.

Uno de estos casos es el de ujo potencial, al cual se llega despreciando los efectosviscosos del uido. D'Alambert (1717-1783) fue una de las primeras personas que aplicó estaformulación. Esto se puede realizar porque en una región alejada de un cuerpo predominaránlos términos convectivos frente a los viscosos. Los ujos que cumplen estas características sonconocidos como irrotacionales.

De esta manera, para realizar las simplicaciones de parte de las ecuaciones 1.1 y 1.2.El hecho de despreciar los esfuerzos viscosos implica que la teoría de ujo potencial debe seraplicada a ujos con valores del número de Reynolds alto, de hecho, lo ideal es aplicarlo paranúmeros de Reynolds innitos.

Como su propio nombre indica, el ujo potencial presenta un campo de velocidadesque deriva de un potencial. Esto es consecuencia directa de que el ujo es irrotacional, porlo que el rotacional del campo de velocidades será nulo:

rot~v =

∣∣∣∣∣∣x y z∂∂x

∂∂y

∂∂z

vx vy 0

∣∣∣∣∣∣ = (−∂vy∂z

,∂vx∂z

,∂vx∂y− ∂vy∂x

) = ~0 (2.1)

Como el estudio se realizará en dos dimensiones el campo de velocidades cumple:

∂vx∂y− ∂vy∂x

= 0

9

CAPÍTULO 2. TEORÍA DEL FLUJO POTENCIAL

El hecho de que el rotacional del campo de velocidades sea cero tiene consecuencias impor-tantes como que el rotacional del gradiente de una función potencial seguirá siendo nulo:

rot∇φ = 0

Por lo tanto se puede denir la función potencial de la velocidad φ, de tal manera que:

−→v = ∇φ

~v = (vx, vy) = (∂∅∂x

,∂∅∂y

)

Por ser irrotacional la circulación a lo largo de cualquier línea cerrada es nula, por lotanto, se podrá aplicar la teoría potencial en las zonas donde no haya vorticidad.

Como se va a resolver para el caso de líquidos, ha de cumplir la ecuación de continuidad1.1:

∂2∅∂x2

+∂2∅∂y2

= 0

o lo que es lo mismo∆∅ = 0 (2.2)

El hecho de que su laplaciano sea nulo se debe a la incompresibilidad del ujo y a queel uido es ideal y rot~v = 0.

Por otro lado, como el uido se considera de densidad constante cumple la ecuación 1.1,o lo que es lo mismo, se puede decir que el ujo es solenoidal. Se puede denir entonces unafunción de corriente tal que su divergencia sea siempre nula. Esta función se llama funciónde corriente ϕ y cumple para el caso plano [11]:

~v = (vx, vy) = (∂ϕ

∂y,−∂ϕ

∂x)

Es sencillo comprobar si se impone que además el ujo sea irrotacional se obtiene:

∂vx∂y− ∂vy∂x

=∂2ϕ

∂y2+∂2ϕ

∂x2= 0

∆ϕ = 0 (2.3)

Será de gran interés el estudio del ujo potencial a través de la función de corrientedado que las líneas de función de corriente constante son líneas de corriente de ujo. Ademáséstas son perpendiculares a las líneas de potencial constante φ.

10 ETSII-UPM


Cabe ahora comentar ahora la relación entre φ y ϕ , lo cual se hará a través de lasiguiente técnica matemática que aparece en las referencias [14] y [11] y que se explica acontinuación.

Se construye una función analítica compleja F (z) tal que su parte imaginaria es lafunción de corriente ϕ y la real es la función potencial φ. Además la variable compleja esz = x+ iy :

F = φ+ iϕ (2.4)

Derivando, por un lado se tiene:

∂F

∂x=dF

dz

dz

dx=dF

dz

∂F

∂y=dF

dz

dz

dy= i

dF

dz

Por otro lado si se deriva en 2.4:

∂F

∂x=∂φ

∂x+ i

∂ϕ

∂x=dF

dz

∂F

∂y=∂φ

∂y+ i

∂ϕ

∂y= i

dF

dz

Por lo tanto, igualando las partes reales e imaginarias de las expresiones anteriores setiene:

∂ϕ

∂y=∂φ

∂x

−∂ϕ∂x

=∂φ

∂y

Como ya se ha comprobado anteriormente, ambas funciones tienen laplaciano nulo y

~v = (vx, vy) = (∂ϕ

∂y,−∂ϕ

∂x) = (

∂φ

∂x,∂φ

∂y) (2.5)

lo cual es congruente con lo que se ha presentado anteriormente. Por la relación matemáticaentre ϕ y φ, ambas funciones son perpendiculares como se observa en la gura 2.1.

Finalmente, se recuerda que la ecuación a resolver en esta parte de la problemática es:

∆ϕ = 0

Por tanto, se ha obtenido una ecuación muy sencilla, en la que la incógnita es la funciónde corriente ϕ, a partir de la cual se conocerá el campo de velocidades. Por otro lado, haciendonulo el término de efectos viscosos en la ecuación 1.2 se obtiene:

ρ~v∇~v = −∇p



Figura 2.1: Líneas de potencial constante y líneas de corriente ortogonales

ρvdv

dx= −dp

dx

Integrando se obtiene un caso particular de la ecuación de Bernouilli:

p

ρ+v2

2= C (2.6)

donde v es el módulo de la velocidad, y C una constante, igual para todo el campo uido.De esta ecuación se obtendrá el campo de presiones. Como ya se ha dicho, esta simplicaciónes válida en zonas alejadas de la geometría a estudiar, ya que en zonas cercanas al cuerpo losresultados experimentales dieren de la teoría de ujo ideal o potencial.

2.2. Condiciones de contorno

La ecuación 2.3 no incluye derivadas temporales, por lo que no serán necesarias condi-ciones de contorno iniciales para su resolución. En cada momento el campo de velocidades seajusta a las condiciones de contorno en ese instante, por lo que la historia previa del uidono inuye en el campo de velocidades posterior del mismo [11]. Esta idea se ve reejada enla imagen 2.2.

En la supercie del sólido la velocidad normal del uido relativa al sólido ha de ser nula,pues físicamente el uido no puede penetrar el sólido. Por lo tanto, el uido solo puede moversetangencialmente a la pared, y ésta constituirá una línea de corriente. Consecuentemente la

12 ETSII-UPM


función de corriente a lo largo de un cuerpo es constante, ya que la supercie del cuerpo esuna línea de corriente como ya se ha comentado. De esta manera, en la geometría a estudiarse impondrá

ϕ|superficie = c

donde c es una constante.

2.3. Paradoja de D'Alambert

Para el caso de ujo irrotacionales la fuerza en el sentido paralelo al ujo FD es nula.Por otro lado, la fuerza en el sentido normal al ujo FL que aparece sobre un cuerpo es [14]

FL = ρΓU0

donde U0 es la velocidad de la corriente incidente y Γ es la circulación. Si la circulaciónes nula, la fuerza en sentido vertical también lo será. Por lo tanto la teoría de ujo potencialpermite llegar a un resultado que, evidentemente, se contradice con la realidad, ya que uncuerpo sumergido en una corriente sí experimenta una fuerza FD distinta de cero. Esto esequivalente a decir que al drag (D) explicado en el Capítulo I es nulo.

En 1752 D'Alambert llegó a la contradicción de que la resistencia producida por uncuerpo inmerso en un uido no puede ser nula, como predice la teoría de ujo potencial. Denuevo, esto signica que el ujo potencial se puede aplicar en regiones alejadas del cuerpo,donde los efectos viscosos son despreciables. En la realidad se observa una resistencia existentepor parte del cuerpo, que se explica debido a la presencia de la capa límite en la región máscercana al cuerpo, como señaló Prandtl en 1904.

Figura 2.2: Flujo ideal



Figura 2.3: Flujo real

14 ETSII-UPM

Capítulo 3

Teoría de la capa límite laminar

En 1904 Prandtl observó que el uido alrededor de un cuerpo se puede dividir en dosregiones. Aquella que queda alrededor del cuerpo en cuestión es la denominada capa límite,en la que los efectos viscosos son predominantes, presentando un gradiente de velocidadesgrande. Teniendo en cuenta las ecuaciones obtenidas 1.1 y 1.2, se procede a la simplicaciónde estas ecuaciones para su posterior tratamiento. La ecuación 1.2 es una ecuación vectorialque se puede descomponer en dos ecuaciones escalares según los ejes X y Y :

ρ(u∂u

∂x+ v

∂u

∂y) = −∂p

∂x+ µ(

∂2u

∂x2+∂2u

∂y2)

ρ(u∂v

∂x+ v

∂v

∂y) = −∂p

∂y+ µ(

∂2v

∂x2+∂2v

∂y2)

3.1. Simplicación de las ecuaciones por análisis dimen-

sional

Si se tiene en cuenta que las ecuaciones de la capa límite se van a aplicar a un uido delongitud característica D en el que el espesor de la capa límite es δ, por análisis dimensionalde los términos se pueden simplicar las ecuaciones. Para ello hay que adimensionalizar todaslas incógnitas, lo cual se hará con las variables de adimensionalización D y U , donde U es lavelocidad de la corriente general. Las ecuaciones anteriores adimensionalizadas, y teniendo encuenta que Re = ρUD

µ= UD

ν. Además es sabido que el espesor de la capa límite es proporcional

a la viscosidad cinemática del uido δ ∼ ν. Por lo tanto, las ecuaciones adimensionalizadasquedan:

u∂u

∂x+ v

∂u

∂y= −∂p

∂x+

1

Re(∂2u

∂x2+∂2u

∂y2)

15

CAPÍTULO 3. TEORÍA DE LA CAPA LÍMITE LAMINAR

u∂v

∂x+ v

∂v

∂y= −∂p

∂y+

1

Re(∂2v

∂x2+∂2v

∂y2)

El análisis dimensional es el siguiente: u ∼ 1, x ∼ 1, v ∼ δ, y ∼ δ, Re ∼ 1δ2. Despre-

ciando los términos pertinentes mediante el estudio de las dimensiones de cada uno se llegaa las ecuaciones de la capa límite:

ρ(u∂u

∂x+ v

∂u

∂y) = −∂p

∂x+ µ

∂2u

∂y2

∂p

∂y= 0

Como la presión a lo largo de la vertical no cambia, ésta será la misma en la zona dela capa límite y en la de ujo potencial, de ahí la validez del acoplamiento. Cabe comentarque estas ecuaciones se pueden aplicar a capas límite laminares, pues si fuera turbulentas elmodelo sería más complejo.

Además, a partir de la ecuación de Bernouilli (2.6) se puede conocer el campo depresiones:

∂p

∂x=dp

dx= −ρU dU

dx

donde U es la velocidad del uido fuera de la capa límite, es decir, de la corriente general.

Por lo general, el espesor de la capa límite aumenta con la longitud del cuerpo en ladirección de la velocidad incidente.

Finalmente, las ecuaciones con las que se va a trabajar serán:

∂u

∂x+∂v

∂y= 0 (3.1)

ρ(u∂u

∂x+ v

∂u

∂y) = −∂p

∂x+ µ

∂2u

∂y2(3.2)

3.2. Condiciones de contorno

Para la resolución de las ecuaciones de la capa límite se impondrán como condicionesde contorno:

u(y = 0) = 0 (3.3)

v(y = 0) = 0 (3.4)

u(y =∞) = U(x) (3.5)

16 ETSII-UPM


u(x0, y) (3.6)

donde U(x) es la velocidad en sentido horizontal de la corriente general y u(x0, y) es elperl de velocidades en la primera sección a estudiar, que será conocido. Las dos primerascondiciones de contorno responden a la condición de no deslizamiento sobre la pared delcuerpo, que se ha explicado en el Capítulo I.

3.3. Desprendimiento de la capa límite

La capa límite no crece indenidamente, llegando a desprenderse en un momento dado.Este fenómeno provoca una disipación de energía que no es deseable en ningún caso. El uidode forma natural avanza, si bien existe un gradiente adverso de presión dp

dx> 0 causado por

lo efectos viscosos que se lo impide:

∂τ

∂y|pared = µ

∂2u

∂y2|pared = ρ(u

∂u

∂x+ v

∂u

∂y)|pared +

∂p

∂x|pared =

∂p

∂x|pared = −ρU dU

dx=dp

dx

donde τ es el esfuerzo cortante. Este gradiente de presiones adverso frena al uido yprovoca el crecimiento de la capa límite. En el punto de desprendimiento se cumple:

∂τ

∂y|pared = 0

En la gura 3.1 puede verse la secuencia que ocurre desde la formación de la capa límitehasta su desprendimiento, donde PI indica el punto de inexión. En el gráco del centro seproduce la separación.

Figura 3.1: Efecto del gradiente de presión adverso



Cabe comentar que desde el momento del desprendimiento de la capa límite en adelante,la teoría de la capa límite que se va a estudiar no es válida, ya que no permite predecir elcomportamiento de ésta de forma real. Habría que resolver las ecuaciones de Navier Stokescompletas.

Figura 3.2: Capa límite adherida y desprendida

3.4. Espesor de la capa límite para la placa plana: ecua-

ción de Karman

Sea el caso de una placa plana, corriente estacionaria y ujo incompresible, la resistenciaal rozamiento viene dada por:

Wr =

ˆ l

0

bτdx

donde τ = µ∂u∂y|y=0 y b es la anchura de la placa. De esta manera, la fuerza de rozamiento

en una sección x de la placa será:

Wr(x) =

ˆ x

0

bτdx

Por otro lado, la fuerza de arrastre D(x), explicada en el Capítulo I, fue obtenida porvon Karman en 1921:

D(x) = ρb

ˆ δ(x)

0

u(U − u)dy = ρbU2θ

donde U es la velocidad de la corriente de la solución potencial, δ(x) el espesor de lacapa límite. Además se tiene el siguiente parámetro:

θ =

ˆ δ

0

u

U(1− u

U)dy (3.7)

18 ETSII-UPM


De esta maneradD(x)

dx= ρbU2 dθ

dx

Igualando ambas expresiones se tiene que:

D(x) =

ˆ x

0

bτdx

dD(x)

dx= bτ

Finalmente, igualando las dos ecuaciones queda:

τ = ρU2 dθ

dx

Von Karman supuso que el perl de velocidades en la capa límite tenía una formaparabólica u(x, y) = U(2y

δ− y2

δ2), e integrando en 3.7 queda θ = 2

15δ y τ = 2µU

δ.

2µU

δ= ρU2 2

15

dδ

dx

Esta ecuación diferencial se integra de 0 a x y de 0 a δ(x):

δ

x=

5,5√Re

Esta ley es la que da el espesor de la capa límite laminar para el caso de una planaplana en función de la longitud x y del número de Reynolds.

3.5. Coeciente de fricción

El coeciente de fricción Cf es la fuerza de rozamiento adimensionalizada, como seexplicó en el Capítulo I En este caso la fuerza a adimensionalizar es la causada por el esfuerzocortante:

Cf =τ

12ρU2

(3.8)

Siguiendo el razonamiento de von Karman se llega a que [14]:

Cf =0,73√Re

Cabe comentar que si la capa límite se desprende, como ya se ha dicho, se cumplirá queτ = µ∂u

∂y|y=0 = 0, o lo que es lo mismo, Cf = 0.



3.6. Solución analítica

En este apartado se estudiar la solución analítica de las ecuaciones de la capa límitepara su posterior aplicación a una placa plana.

Sea un ujo que pasa a través de una cuña con ángulo πβ, retomando las ecuaciones dela capa límite y sus condiciones de contorno 3.1, 3.2, 3.3, 3.4, 3.5 y 3.6, se realiza el siguientecambio de variable para adimensionalizar y, el cual se sugiere en la bibliografía [12]:

η = y

√m+ 1

2

U(x)ρ

xµ(3.9)

πβu1

Figura 3.3: Cuña de ángulo πβ

En este caso, según sugería la bibliografía consultada [12] se ha tenido en cuenta que lacorriente de ujo potencial, exterior a la capa límite, lleva una velocidad U(x) = u1x

m dondeu1 es una constante y m es una parámetro relacionado con el ángulo de la cuña:

m =β

2− β (3.10)

De esta manera, el sistema a resolver quedará:

f ′′′ + ff ′′ + β(1− f ′2) = 0 (3.11)

f(η = 0) = 0

f ′(η = 0) = 0

f ′(η =∞) = 1

las velocidades horizontal y vertical en la capa límite quedan:

u = U(x)f ′(η)

v = −√m+ 1

2

µ

ρ

U(x)

x(f(η) +

m− 1

m+ 1ηf ′(η))

20 ETSII-UPM


Resolviendo las ecuaciones para distintos valores de m se obtienen las grácas queaparecen en la gura 3.4.

Figura 3.4: Curvas de Falker-Skan. Fuente [12]

3.6.1. Caso: placa plana

Para este caso particular, la solución a la ecuación 3.11 es la solución de Blasius (1908):

η =y

δ= y

√U(x)ρ

xµ= y

√u1ρ

xµ(3.12)

u = U(x)f ′(η) = u1f′(η)

v = −√µ

ρ

u1x

(f − ηf ′)

Realizando estos cambios de variable en las ecuaciones de la capa límite se obtiene laecuación resultante:

f ′′′ +1

2ff ′′ = 0 (3.13)

Se ha resuelto esta ecuación diferencial con Matlab para poder compararla posteriormentecon la solución obtenida a través de la resolución de las ecuaciones de la capa límite laminarcon condiciones de contorno de ujo potencial para el caso de una placa plana.

Por lo tanto, si u es la velocidad dentro de la capa límite y U lo es fuera, para unaaltura suciente de la capa límite se tendrá que u

U' 1. Según la bibliografía consultada esto

ocurre para η ' 5, por lo que el espesor de la capa límite será:



δ = 5

√µx

ρu1=

5x√Re

(3.14)

Es de especial interés obtener el coeciente de rozamiento Cf , que se recuerda que es la fuerzade rozamiento adimensionalizada:

Cf =τ

12ρU2

τ = µ∂u

∂y|y=0 = µU

∂f ′

∂η

∂η

∂y|y=0 = µU

√Uρ

xµf ′′(0)

De la bibliografía consultada se obtiene que f ′′(0) = 0,332 [12], dato que también seobtendrá posteriormente en Matlab. De esta manera el coeciente de fricción será:

Cf =τ

12ρU2

=0,332µU

√Uρxµ

12ρU2

=0,664√Re

(3.15)

Como se puede observar, es muy aproximado al que predijo Karman Cf = 0,73√Re.

3.6.2. Solución pseudoanalítica

En este caso se va a resolver por métodos numéricos la ecuación que da la soluciónanalítica para una placa plana (3.13). Las condiciones de contorno son f(η = 0) = 0, f ′(η =0) = 0 y f ′(η =∞) = 1.

Por lo tanto, se ha obtenido la solución analítica con un solver de Matlab para la resolu-ción de ecuaciones lineales en derivadas parciales. Este solver es ode45 , basando en el métodode Runge-Kutta para la resolución de ecuaciones diferenciales, que se sale del temario abar-cado en este proyecto. Si bien, este método sirve para aproximar las ecuaciones diferencialesa las de soluciones de ecuaciones diferenciales ordinarias. De esta manera, se crea la función@placa, que pone de forma explícita la ecuación diferencial. Los argumentos del solver ode45son la función @placa, el intervalo de η donde se va a resolver la ecuación diferencial, y lascondiciones de contorno expuestas más arriba. El solver ode45 realmente no permite obtenerla solución analítica ya que resuelve la ecuación diferencial por métodos numéricos, si bienla solución es muy aproximada. La llamaremos por ello solución pseudoanalítica.

Como ya se ha comentado, para η ' 5 se alcanzaba la altura máxima de la capalímite, resolviéndose, por tanto, todo el campo uido de la capa límite. Sin embargo, paraasegurarnos de que se resuelve la capa límite completa en altura, el vector que recorre éstaen la dirección vertical irá hasta 7. No hay ningún problema en resolver más allá de la capalímite, ya que el perl de velocidades de Blasius, que puede verse en la gura 3.7 llega a

22 ETSII-UPM


velocidad 1 a cierta altura próxima a 5, y más arriba de ese punto la velocidad es constantee igual a 1, pues fuera de la capa límite se tiene ujo potencial.

Posteriormente se comprueba que uU' 1 en η = 7. Por lo tanto queda justicada la

validez de resolver para un vector η de 0 a 7, con intervalos de 0,01, valor que también se hadeterminado experimentalmente, ya que la solución obtenida no varía para mayores renados,es decir, para intervalos de η menores de 0,01.

El tercer argumento del solver ode45 es un vector que contiene las condiciones decontorno de f , f ′ y f ′′ en η = 0. Las condiciones de contorno sobre los dos primeros sontriviales, al venir dadas por el planteamiento inicial del problema, si bien la imposición dela última condición de contorno no es posible analíticamente dado que η =∞ se reere al ηdonde termina la capa límite, además de que se necesita una condición de contorno sobre f ′′.Para resolver este problema se ha recurrido el método del disparo.

El método del disparo impone a través de un bucle distintas condiciones de contornopara f ′′(0) y resuelve la ecuación diferencial. Finalmente se prueba para cuál de estos valoresse ha cumplido que f ′(η =∞) = 1. La gura 3.5 muestra el resultado obtenido, que es:

f ′′(η = 0) = 0,332

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Condición de contorno

0

0.5

1

1.5

2

2.5

df dη| ∞

Figura 3.5: Condición de contorno sobre f ′′

Con esta condición de contorno ya conocida se procede a la obtención de la soluciónpseudoanalítica para obtener el conocido perl de velocidades en la capa límite laminar para



el caso de una placa plana. También se obtendrá la curva de f ′ en función de η, coincidentecon la mostrada en la gura 3.4 para m = 0. Los resultados obtenidos, por tanto, aparecenen la gura 3.6. La solución que aparece en la gura 3.7 fue la obtenida por Blasius, y queaparece en la bibliografía consultada [12].

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

η

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

df dη

Curva teórica de la distribución de velocidades en la capa límite laminar

Figura 3.6: Curva obtenida, igual a las de Falker-Skan

24 ETSII-UPM


0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

U

0

1

2

3

4

5

6

7

Altu

ra

Perfil de velocidades en la capa límite para una placa plana

Figura 3.7: Perl de velocidades obtenido: perl de Blasius



26 ETSII-UPM

Capítulo 4

Modelo matemático

El presente capítulo tiene como objetivo la puesta en práctica de la metodología adesarrollar. Para ello se explicará de manera general para aplicarlo en los siguientes capítulosa los casos de la placa plana y del cilindro.

4.1. Introducción a los métodos numéricos

Los métodos numéricos son una herramienta matemática que permite simplicar laresolución de problemas, obteniendo soluciones aproximadas a las analíticas. A menudo seutilizan en la resolución de problemas que no tienen solución analítica o es muy compleja suobtención. Se hará uso del análisis numérico en este proyecto con el n de resolver el campode velocidades tanto en la zona de la capa límite como en la de ujo potencial.

Cabe comentar que la resolución de ecuaciones diferenciales por métodos numéricosutilizada requiere una primera adimensionalización de las mismas. Posteriormente se continúaaproximando las ecuaciones por el método de elementos nitos o por el métodos de diferenciasnitas.

El método de elementos nitos es el que utiliza el programa FreeFem++ [1] para resolverla ecuación del ujo potencial, mientras que el de diferencias nitas es el que se ha aplicadopara resolver las ecuaciones de la capa límite laminar en Matlab [2].

4.2. Adimensionalización de las variables

Como ya se ha comentado, el número de Reynolds es:

Re =ρU∞D

µ

27

CAPÍTULO 4. MODELO MATEMÁTICO

A continuación se procede a la adimensionalización de las variables teniendo en cuentaque las variables para adimensionalizar son ρ, U∞ y D, donde U∞ es la velocidad de lacorriente, D es una longitud característica del problema.

Las nuevas variables adimensionalizadas son las que llevan comilla:

p′ =p

ρU2∞

~v′ =~v

U∞

x′ =x

D

∂′

∂′x′=

∂

∂xD

Tanto las variables utilizadas en FreeFem++ como en Matlab son adimensionales, porlo que las ecuaciones introducidas en los programas también lo son.

4.3. Datos del problema

Los datos comunes a los problemas de ujo potencial y capa límite son la geometría, elnúmero de Reynolds y la velocidad de la corriente en el innito U∞.

Cada problema tendrá una geometría determinada, y el mallado de cada problema segenerará en base a ésta. Como la variable de adimensionalización es U∞, la corriente en elinnito adimensionalizada será U∞ = 1 para todas las geometrías y para los dos problemasde vericación del método desarrollado.

En cuanto al número de Reynolds, se recuerda que lo ideal es tener en teoría de ujopotencial el Reynolds innito, lo cual no va a ser posible; y para las ecuaciones de la capalímite un número de Reynolds laminar. La resolución del campo potencial no requiere laintroducción de un número de Reynolds, por lo que el Re elegido será para la capa límite.Como se verá más adelante, las ecuaciones de la capa límite se adimensionalizan de tal maneraque no dependen del número de Reynolds. Así, lo único que hay que tener en cuenta es que losresultados obtenidos en la parte de ujo potencial son válidos para un número de Reynoldselevado.

Las condiciones de contorno se obtienen de estos datos para ambos problemas.

28 ETSII-UPM


4.4. Campo de velocidades según la teoría del ujo po-

tencial.

Se recuerda que la ecuación a resolver es

∆ϕ = 0 (4.1)

siendo las condiciones de contorno

ϕ|superficie = c

donde c es una constante y ϕ función de corriente; y otras condiciones sobre las velo-cidades teniendo en cuenta que

(vx, vy) = (∂ϕ

∂y,−∂ϕ

∂x)

Para cada geometría, en concreto, se deducirán las condiciones de contorno sobre ϕ.Además será dato la velocidad de la corriente aguas arriba U∞.

la resolución de la ecuación 4.1 se hace por elementos nitos, con el programa Free-Fem++ [1]. Previamente será necesario acondicionar la ecuación 4.1 para poder aplicar lateoría de métodos numéricos en este problema

Estos elementos nitos estarán en el espacio de elementos nitos Vh. Este espacio es unespacio de dimensión nita que contiene a las funciones ϕ y g, que son las funciones base,y que está denido en todos los elementos triangulares del mallado. En el espacio Vh lasfunciones base g son polinomios de grado m en cada elemento triangular, y satisfacen que enla frontera de Dirichlet siempre toman valor 0. En cambio, las funciones base ϕ, que son lasolución del problema, tomarán un valor determinado en cada frontera.

De esta manera, la ecuación 4.1 habrá que multiplicarla por las funciones base de g eintegrar en el dominio Ω: ˆ

Ω

4ϕ g dΩ = 0

Integrando por partes queda:ˆΩ

∇ϕ∇g dΩ −ˆ∂Ω

∂ϕ

∂n·g dΓ = 0

Sobre la frontera ∂Ω la función base g es nula como ya se ha dicho, por lo que se anulala segunda integral en toda la frontera. Hacer nula las funciones base en la frontera es una



condición de contorno que se impone para reducir el número de ecuaciones a resolver, esdecir, se dan por supuesta las funciones base g en los nodos de la frontera. Por lo tanto laecuación a resolver en FreeFem++ queda:

ˆΩ

∂ϕ

∂x

∂g

∂x+∂ϕ

∂y

∂g

∂ydΩ = 0

Finalmente, las condiciones de contorno a imponer serán de tipo Dirichlet, es decir, ϕserá conocida en todas las fronteras.

4.4.1. Mallado

La resolución en FreeFem++ requiere la denición de la geometría y el mallado. Paraello se dene la geometría del problema, que se explicará más adelante para cada uno de ellos,y el número de elementos que se quieren en cada frontera. A partir de ahí es el programael que crea el mallado interno de la geometría, con el número de elementos que encajen enesa disposición. Además existen funciones de FreeFem++ que realizan un proceso iterativode resolución del campo de ujo potencial y adaptación del mallado a la misma, pudiendointroducir el máximo tamaño de celda del mallado o el error máximo permitido. Esta funciónse llama adaptmesh. A modo de ejemplo se muestra en la gura 4.1 un ejemplo de malladoalrededor de un perl NACA.

Figura 4.1: Ejemplo de mallado en FreeFem++

4.4.2. Discretización y datos

FreeFem++ resuelve internamente la ecuación por elementos nitos, por lo que nohace falta las ecuaciones sean discretizadas. Lo único que hay que tener en cuenta es que lasvariables dato han de estar adimensionalizadas, según se ha explicado en la sección 4.2.

30 ETSII-UPM


Como ya se ha comentado, los datos de entrada a este problema será la velocidad enel innito, que es 1 adimensionalizada, así como las condiciones de contorno de Dirichlet encada frontera de la geometría.

4.4.3. Resolución

Por lo tanto, FreeFem++ obtiene la función de corriente ϕ en cada nodo del mallado, yderivando se obtiene el campo de velocidades (vx, vy) adimensionalizadas. De igual manera,se puede obtener el campo de presiones p a través de la ecuación de Bernouilli (2.6), y sugradiente ∂p

∂x.

Si se considera la presión nula en corriente arriba, lo cual se ha considerado en todoslos problemas resueltos, donde la velocidad es U∞ entonces la ecuación 2.6 queda:

U2∞2

= C

Introduciendo el valor de la constante C de nuevo en la ecuación 2.6 para un puntogenérico se obtiene:

p =ρ

2(U2∞ − v2) =

ρ

2(U2∞ − v2x − v2y)

Si bien el campo de presiones que se obtiene en FreeFem++ es adimensional, de talmanera que:

p =ρ

2(U2∞ − v2)

ρU2∞p′ =

ρ

2(U2∞ − U2

∞v′2)

p′ = 0,5(1− v′2)

Retomando la notación sin comillas, se tiene la presión adimensionalizada, que será lacalculada en FreeFem++ : p = 0,5(1− v2).

4.5. Campo de velocidades en la capa límite laminar

Se recuerda que las ecuaciones de la capa límite laminar y sus condiciones de contornoson:

∂u

∂x+∂v

∂y= 0



u∂u

∂x+ v

∂u

∂y= −1

ρ

∂p

∂x+µ

ρ

∂2u

∂y2

u(y = 0) = 0

v(y = 0) = 0

u(y =∞) = U(x)

u(x0, y)

donde u(x0, y) es la velocidad de la corriente al inicio de la geometría. El problema seresuelve en Matlab.

Es necesario adimensionalizar las ecuaciones de la capa límite, si bien por simplicidadse va a adimensionalizar la ecuación vectorial de Navier Stokes para deducir las ecuacionesde la capa límite laminar adimensionalizadas.

ρ~v∇~v = −∇p+ µ4~v

Por lo tanto, teniendo en cuenta la adimensionalización de las variables explicadas enla sección 4.2:

ρU∞~v′U∞D∇′~v′ = −ρU

2∞D∇′p′ + ρU∞D

Re

U∞D24′~v′

Simplicando queda

~v′∇~′v′ = −∇′p′ + 1

Re4~′v′

Es importante tener en cuenta que, para comodidad del lector, las variables adimen-sionalizadas se van a tratar por su nombre anterior sin adimensionalizar, es decir, se quitanlas comillas en las variables adimensionales. Para evitar confusiones más adelante, todas lasvariables que aparezcan serán adimensionales a no ser que se comente lo contrario. Por lotanto, retomando la notación anterior, la ecuación de la capa límite adimensionalizada queda:

∂u

∂x+∂v

∂y= 0

u∂u

∂x+ v

∂u

∂y= −∂p

∂x+

1

Re

∂2u

∂y2(4.2)

32 ETSII-UPM


4.5.1. Mallado

Será necesaria la discretización de las ecuaciones de la capa límite, para lo cual hayque tener en cuenta el mallado. Para ello será necesario conocer la geometría del problemaa resolver, si bien Matlab [2] hace uso de este mallado para el planteamiento del sistema deecuaciones con incógnitas que se resuelve, como se explica más adelante. De esta manera, elmallado de Matlab es cticio, pues es el que se impone para construir el sistema de ecuaciones.Este mallado es regular como el que se muestra en la gura 4.2, y es el que permite el cálculodel sistema de ecuaciones en Matlab.

Por lo tanto, cualquier elemento diferencial de la supercie de un cuerpo se aproximaa una placa plana de longitud diferencial. El mallado a utilizar será tal que el eje X tendráNx nodos, indicados por el subíndice n, separados una distancia de 4x mientras que el ejeY tendrá N nodos, separados una distancia 4y, y numerados con el subíndice j .

x

y

n n+ 1n− 1

j + 1j

j − 1

Up

Figura 4.2: Mallado para los elementos n− 1, n y n+ 1

4.5.2. Discretización y datos

En primer lugar, cabe comentar que la discretización se puede realizar por diferenciasnitas o por elementos nitos. Como ya se ha comentado, se ha elegido el método de dife-rencias nitas, discretizando según cada eje. Las derivadas parciales se han aproximado pordiferencias nitas donde el error cometido es de orden O(4x), como se muestra a continuaciónteniendo en cuenta la nomenclatura de la gura 4.2.



Discretización según el eje X:

∂u

∂x'unj − un−1j

4x

Discretización según el eje Y :

∂v

∂y' vnj − vnj−1

4y∂u

∂y' unj − unj−1

4y∂2u

∂y2' unj+1 − 2unj + unj−1

4y2

Finalmente, las ecuaciones a resolver quedan:

unj − un−1j

4x +vnj − vnj−14y = 0 (4.3)

un−1j

unj − un−1j

4x + vn−1j

unj − unj−14y = −∂p

∂x|n +

1

Re

unj+1 − 2unj + unj−14y2 (4.4)

donde j = 1, 2, 3, ...N y n = 1, 2, 3, ...Nx.

Las condiciones de contorno a imponer serán las presentadas al inicio de esta sección,que transformadas quedan:

un1 = 0

vn1 = 0

unN = Unp

u1(y) = u(x0, y)

pn = pnp

para todo n. Las dos primeras condiciones de contorno imponen la condición de nodeslizamiento. La tercera impone que la velocidad en sentido horizontal fuera de la capalímite (nodo (n,N)) sea la de la solución potencial Un

p . Además, el perl de velocidades paran = 1 debe ser conocido, lo cual queda reejado en la cuarta condición de contorno. La quintacondición de contorno es la presión pnp obtenida del campo de ujo potencial para cada puntodel mallado n.

34 ETSII-UPM


La tercera condición de contorno se cambia, para algunos casos, por:

unN = unN−1 (4.5)

al dar mejores resultados que la anterior. Esto se ha realizado habiendo probado laanterior a esta y comprobando que ambas velocidades son prácticamente iguales a Un

p . Estacondición de contorno viene de evaluar la ecuación 4.2 en el nodo N , de tal manera que seobtiene:

u(y = δ)∂u

∂x|y=δ + v(y = δ)

∂u

∂y|y=δ = −∂p

∂x|y=δ +

1

Re

∂2u

∂y2|y=δ

Como u(y = δ) = Up y fuera de la capa límite la velocidad es horizontal, entonces∂2Up

∂y2|y=δ = 0. Además, en una aproximación al ujo de Bernouilli (ujo ideal) es sabido que

− ∂p∂x|y=δ = Up

∂Up

∂x. De esta manera se obtiene que ∂Up

∂y|y=δ = 0, por lo que queda justicada la

condición de contorno 4.5.

Por lo tanto, el acoplamiento entre los problema de ujo potencial y capa límite serealiza a través de − ∂p

∂x|n, que será conocido en todos los puntos del mallado. En el Capítulo

III se dedujo que ∂p∂y

= 0, por lo que la presión a lo largo de la vertical es constante e igualpara el campo de ujo potencial y la zona de la capa límite. De esta manera, la resoluciónen FreeFem++ permite obtener un vector de presiones p, que se importará a Matlab, siendoconocido − ∂p

∂x|n en cada punto.

4.5.3. Resolución

A continuación se procede a explicar la resolución del problema a partir de las ecuacio-nes, los datos, y las condiciones de contorno. El objetivo es conocer el campo de velocidadesen la capa límite laminar, entre otras cosas, para a partir de él, obtener los perles a lo largodel cuerpo y parámetros interesantes como el coeciente de fricción Cf .

En Matlab se ha realizado un programa que obtiene las matrices U y V , donde U es lavelocidad horizontal en la capa límite y V la vertical. Estas matrices contienen los valores delas velocidades u y v en cada nodo del mallado. Se resolverán de forma unidimensional, de talmanera que para cada n se obtendrán Un y V n, y una vez conocidos el programa resolverála columna de nodos n+ 1.

Para mayor comodidad del lector a la hora de entender el programa se van a ordenarlas ecuaciones en función de las incógnitas que contienen para la posterior construcción delas matrices que conforman el sistema de ecuaciones.

En función de la geometría del problema se construyen los vectores x e y, que recorrenel mallado con intervalos de 4x y 4y. Los valores de estos incrementos se han elegido trasprobar diferentes valores para cada geometría. Como cualquier estudio de métodos numéricos,



es interesante la comparación de resultados para distintos pasos del mallado, si bien en esteproyecto la elección del mallado adecuado se realiza por análisis de los perles de velocidadesy vericación de los mismos. Por ejemplo, para el caso de la placa plana, se tiene que obtenerel perl de velocidades de Blasius, por lo que se probarán distintos valores de 4x y 4y. Paratodas las geometrías resueltas estarán comprendidos entre 0,01 y 0,001.

De esta manera, retomando las ecuaciones 4.3 y 4.4 quedan:

(−vn−1j

4y −1

Re4y2 )unj−1 + (un−1j

4x +vn−1j

4y +2

Re4y2 )unj + (− 1

Re4y2 )unj+1 = −(un−1j )2

4x − ∂p

∂x|n

(4.6)

− 1

4yvnj−1 +

1

4yvnj = −

unj − un−1j

4x (4.7)

Como se puede observar, la primera ecuación tiene como incógnitas unj−1, unj y unj+1.

Todas las variables evaluadas en n − 1 son conocidas ya que vendrán de la resolución de lacolumna de nodos anterior. Una vez resuelta la ecuación primera se procede a la resoluciónde la segunda donde las incógnitas son vnj−1 y v

nj .

Es sabido que el espesor de la capa límite δ es proporcional al número de Reynolds dela siguiente manera:

δα1√Re

Por lo tanto, la altura y también lo será. Si se realiza el cambio de variable y′ =√Rey,

siendo y′ la nueva variable y, se tienen ecuaciones más simplicadas. Lo que se consigue coneste cambio de variable es que la nueva variable y′ sea del orden de 1, además de conseguirque las ecuaciones de la capa límite no dependan del número de Reynolds. Además se realizael cambio de variable vnj ´ = vnj

√Re.

Una vez realizados estos cambios de variable, por simplicidad, se asigna la siguientenotación:

an = −vn−1j

4y −1

4y2

bn =un−1j

4x +vn−1j

4y +2

4y2

c = − 1

4y2

fn = −(un−1j )2

4x − ∂p

∂x|n

36 ETSII-UPM


ln = −unj − un−1j

4x (4.8)

Por tanto, el sistema de ecuaciones 4.6 en forma matricial queda de la siguiente manera:

A =

b1 c 0 · · · 0 0 0 0a2 b2 c · · · 0 0 0 0...

......

. . ....

......

...0 0 0 · · · bN−2 c 0 00 0 0 · · · aN−1 bN−1 c 00 0 0 · · · 0 aN bN c

U =

un1un2...

unN−2unN−1unN

=

0un2...

unN−2unNunN

B =

f1f2...

fN−2fN−1fN

AU = B

Una vez resuelto éste se resuelve 4.7, que en forma matricial es:

K =

14y 0 0 · · · 0 0 0

− 14y

14y 0 · · · 0 0 0

0 − 14y

14y · · · 0 0 0

......

.... . .

......

...0 0 0 · · · 1

4y 0 0

0 0 0 · · · − 14y

14y 0

0 0 0 · · · 0 − 14y

14y



V =

vn1vn2vn3...

vnN−2vnN−1vnN

=

0vn2vn3...

vnN−2vnN−1vnN

L =

l1l2l3...

lN−2lN−1lN

KV = L

Los vectores U y V dan el campo de velocidades en todo el mallado.

4.5.4. Coeciente de fricción

Con los perles de velocidades se puede calcular el valor del coeciente de fricción (3.8)en cualquier punto de la supercie del cuerpo a analizar. Se recuerda que :

Cf =τ

12ρU2∞

=2µ

ρU2∞

∂u

∂y|y=0

Teniendo en cuenta las variables adimensionalizadas y aproximando la derivada pordiferencias nitas queda:

Cf =2µ

ρU2∞U∞

√Re

D

Un2 − Un

1

4y =2√Re

Un2 − Un

1

4y (4.9)

donde D es la longitud característica del problema.

38 ETSII-UPM

Capítulo 5

Flujo alrededor de una placa plana

A continuación se procede a la resolución completa del problema de ujo alrededor deuna placa plana. Este caso fue resuelto analíticamente por Blasius según se ha visto en elCapítulo III. Es el caso de geometría más sencillo, y además se dispone se su solución analítica,lo que lo hace idóneo para compararlo con los resultados obtenidos del modelo matemáticodesarrollado. Por último, se ha introducido uno de los métodos de control de la capa límitemás frecuentes, que es la succión, y que se explicará en el apartado correspondiente.

Para una placa plana, el número de Reynolds de transición entre régimen laminar yturbulento es 105. De esta manera, se ha elegido Re = 104 por ser un número de Reynoldselevado pero aún en régimen laminar. En cuanto a la geometría, la longitud de la placa Lserá la dimensión característica del problema, por lo que la longitud adimensionalizada es 1.La velocidad en el innito, como ya se ha dicho, será U∞ = 1.

Re =ρU∞L

µ

Se espera que la capa límite no se desprenda al no existir un gradiente adverso depresiones, por lo que el coeciente de fricción será mayor que 0 y no nulo para todo eldominio.

Se adelanta además, la trivialidad de la solución del ujo potencial, por lo que se podránresolver las ecuaciones de la capa límite deduciendo las condiciones de contorno procedentesde la resolución del ujo potencial a partir de las características del ujo alrededor de estageometría. Sin embargo, se ha realizado para mostrar cómo funciona esta metodología delforma genérica.

5.1. Flujo potencial

La solución del campo de velocidades de la teoría del ujo potencial es trivial como sedemuestra a continuación.

39

CAPÍTULO 5. FLUJO ALREDEDOR DE UNA PLACA PLANA

Las condiciones de contorno son tales que en la placa la función de corriente es cons-tante, e igual a 0, según se ha elegido. En la frontera superior también se ha tomado comoconstante al estar muy alejada de la placa plana. Esto es posible porque la geometría essucientemente grande como para poder despreciarse los efectos viscosos en todo el dominioy aplicar ujo potencial. Como las líneas de corriente en el campo uido derivado poten-cial son paralelas a la placa plana, la frontera superior es también una línea de corriente, yconsecuentemente la función de corriente ϕ es constante. En este caso, a la entrada se haconsiderado un perl recto y uniforme de velocidad, cuyo valor adimensionalizada es 1 comoya sabíamos:

(vx, vy) = (∂ϕ

∂y,−∂ϕ

∂x) = (1, 0)

Integrando se obtiene una función de corriente que será:ϕ = y.

Por lo tanto, en la frontera superior será la altura del mallado. De igual manera que ala entrada, a la salida también se tiene ϕ = y, pues el perl de velocidades a la salida es elmismo que a la entrada, como el que aparece en la gura 5.1.

Figura 5.1: Perl de entrada a la placa plana

Las dimensiones de este problema no inuyen en su resolución, si bien es necesariodenir la geometría del mismo, por lo que se ha elegido un altura y una longitud determina-das. En este caso las dimensiones de la geometría son irrelevantes porque, se recuerda, queFreeFem++ resuelve la ecuación 4.1, que está adimensionalizada.

Las condiciones de contorno quedan como las de la gura 5.2.

40 ETSII-UPM


C1

C2

C3

C4

ϕ = y

ϕ = 1

ϕ = y

ϕ = 0

Figura 5.2: Condiciones de contorno para ujo potencial en la placa plana

Resolviendo el problema se obtienen las funciones de corriente a lo largo de la geometría,que serán paralelas a la placa plana.

Como se puede observar, la solución es trivial, de tal manera que las funciones decorriente son paralelas al sentido del ujo, lo cual ya era sabido. En todo el campo uido lavelocidad es idénticamente igual a 1 en el sentido horizontal, independientemente del puntoen que se evalúe.

El programa FreeFem++ da un gradiente de presiones del orden de 10−14, que se puedeaproximar a 0, lo cual era de esperar dado que

p =1

2(1− v2) =

ρ

2(1− v2x − v2y) =

1

2(1− 12 − 02) = 0

Por lo tanto, este vector de presiones es nulo de cara a la resolución de la capa límite conMatlab, ya que importar un vector con valores del orden de 10−14 solo dará error numérico.De esta manera, se justica la poca importancia que tiene la resolución del ujo potencial, yconsecuentemente las dimensiones de la placa. Sin embargo, este es un caso aislado, pues lasolución de este apartado es trivial como ya se ha dicho.

5.2. Capa límite

Para esta geometría se ha elegido 4x = 10−4 y 4y = 10−2, tras haber probado condiferentes valores y ver que para estos intervalos los errores numéricos por la discretizaciónen X e Y se hacen despreciables.

Para el caso de la placa plana, la condiciones de contorno a imponer son:

u1(y) = u(x0, y) = 1



pn = pnp = 0

unN = Unp = 1

para todos los nodos n, donde Unp es la velocidad del ujo potencial en el nodo n. Por

lo tanto, el primer perl de velocidades se ha tomado como un perl recto igual que el quese ha introducido como perl inicial para la resolución del ujo potencial.

El mallado ya se ha explicado en la sección 4.5.1. De esta manera, se recuerda que seresolverá para cada nodo de la placa n, una la vertical de nodos j.

Como ya se ha comentado el vector x que recorre la placa en sentido horizontal iráhasta 1, mientras que el vector y irá hasta 7, según lo expuesto en la sección 3.6.2. De estamanera, η (3.12) e y son equivalentes como se demuestra a continuación:

η =y

δ= y

√U(x)ρ

xµ= y

√u1ρ

xµ

Utilizando las variables adimensionales queda:

η = y′√

u1ρ

x′Lµ

L√Re

=y′√x′

donde estas ya son las variables adimensionalizadas. Por lo tanto, si ηmax = 7, para elvalor mayor espesor de la capa límite, que se dará en x = 1, se tiene y = 7.

En primer lugar, para el intervalo horizontal se ha tomado 4x = 10−4, por ser unintervalo adecuado para el estudio de los cambios de la capa límite. El vector x que recorre laplaca plana irá de 0 hasta la longitud de la placa adimensionalizada L

L= 1 Consecuentemente,

el valor del número de nodos en el sentido horizontal será:

Nx =1

4x = 10001

N =7

4y + 1 = 701

Se resuelve el problema y se obtienen los perles de velocidades de la gura 5.3, dondese han representado en cuatro puntos de la placa.

Como se puede observar las velocidades horizontales U a lo largo de la placa crecen enaltura, y valen 1 a partir de una determinada altura. Se distingue el crecimiento de la capalímite conforme avanza en posiciones el ujo. El primer valor representado es el tomado parael primer nodo, e igual a 1 para toda la altura.

42 ETSII-UPM


0 0.5 1

u

0

1

2

3

4

5

6

7

Altu

ra

x=0.25

0 0.5 1

u

x=0.5

0 0.5 1

u

x=0.75

0 0.5 1

u

x=1

Figura 5.3: Perles de velocidades a lo largo de la placa

5.2.1. Comparación con la solución pseudoanalítica

Se recuerda que en el Capítulo III se resolvió la ecuación diferencial que da la soluciónanalítica al problema de la placa plana, obteniéndose una solución pseudoanalítica.

Para vericar la validez de los resultados obtenidos en Matlab se han superpuesto lasgrácas obtenidas por ambos métodos, como se muestra en la gura 5.4.

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

η

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

df dη

Solución pseudoanalíticaSolución numérica

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

U

0

1

2

3

4

5

6

7

Altu

ra


Figura 5.4: Comparación con las solución pseudoanalítica

La solución numérica que se muestra en estas grácas es la correspondiente al últimonodo, donde x = 1, ya que la equivalencia η = y se da en ese punto. Además u = U∞f

′(η),



por lo que se puede establecer fácilmente la relación entre ellos. Como se puede ver en elcódigo en los Anexos, la solución u numérica se ha pasado a f ′ y la solución pseudoanalíticaf ′ se ha pasado a u. Así ha sido posible construir estas grácas.

5.2.2. Coeciente de fricción

En este apartado se va a comparar el coeciente de fricción Cf teórico con el que seha obtenido por métodos numéricos, así como el cálculo del coeciente de fricción medio a lolargo de la placa plana.

Teniendo en cuenta la ecuación 4.9 y conociendo las velocidades horizontales en todoslos puntos del mallado del problema de la capa límite, es posible calcular el coeciente defricción.

Por otro lado, teniendo en cuenta las variables adimensionalizadas en Matlab el teórico(3.15), queda:

Cf (x) =0,332µU

√Uρxµ

12ρU2

=0,332µU

√Uρx′µL

12ρU2

=0,664√Re√x′

(5.1)

En la gura 5.5 se han representado ambos coecientes para su comparación, que seaprecia mejor en ejes logarítmicos, y se ha representado a partir de una distancia de 0,01desde el inicio de la placa.

10-2 10-1 100

Longitud adimensionalizada de la placa

10-2

10-1

100

logC

f(x)

Coeficiente de fricción en ejes logarítmicos


Figura 5.5: Coeciente de fricción en ejes logarítmicos

44 ETSII-UPM


Para los puntos situados a una distancia menos de 0.01 la solución obtenida se desvíamás de la teórica. Esto es porque el perl inicial es recto y el que corresponde a la soluciónde Blasius es parabólico. La transformación del primero en el segundo conlleva una serie deerrores que se ven reejados en el coeciente de fricción. podría contemplarse la opción derenar más en mallado en los primeros nodos de la placa plana, si bien, el método utilizadono permite un renamiento local. Es por ello que se ha tomado 4x sucientemente pequeñocomo para que la solución sea válida como se acaba de ver en la gura 5.4. Además secomprueba que el coeciente de fricción es mayor que 0 en todos los puntos de la placa porexistir su logaritmo. Esto ya se había adelantado al inicio de este capítulo y es indicativo deque no ha habido desprendimiento de la capa límite.

Además se ha calculado el error cometido con respecto al teórico (5.1), como muestrala gura 5.6.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-100

-80

-60

-40

-20

0

20

40Error entre coeficientes de fricción teórico y numérico en %

Figura 5.6: Error cometido en Cf

El error cometido es elevado al principio porque en los primeros puntos de la placa planael perl de velocidades pasa de valer 1 a valer 0 en la pared del cuerpo, como se ha explicadounas líneas más arriba. De hecho las ecuaciones de la capa límite dejarían de ser válidas alcomienzo de la placa plana, pues en esta zona el espesor de la capa límite no es mucho menorque la longitud en esta zona. Este cambio de la geometría en el perl de velocidades se da enel primer intervalo 4x, y es por lo que el error presenta un pico en el primer nodo. Una vezestabilizado este cambio brusco se aprecia que el error cometido es muy bajo.

En cuanto al coeciente de fricción medio, el que se obtiene por métodos numéricos es:

Cf =

´ L0τ dx

12ρU2∞L

=2µ

ρU2∞L

ˆ L

0

∂u

∂y|y=0dx



Teniendo en cuenta las variables adimensionales se tiene:

Cf =2µ

ρU2∞L

U∞

ˆ 1

0

∂u

∂y|y=0

√Re

LLdx =

2√Re

ˆ 1

0

∂u

∂y|y=0dx

donde la integral se ha aproximado por reglas de cuadratura numérica. En este caso seha utilizado la regla del trapecio:

ˆ b

a

f(x)dx ' (b− a)f(a) + f(b)

2(5.2)

y en el caso de cuadratura compuesta:

ˆ 1

0

f(x)dx 'Nx∑k=2

4x2

(f(k) + f(k − 1))

donde f(i) es la derivada ∂u∂y|y=0 evaluada por métodos numéricos:

f(i) ' U i2 − U i

1

4y

Realmente se ha calculado Cf√Re por no depender la resolución de las ecuaciones del

número de Reynolds, y el resultado obtenido ha sido el siguiente:

Cf =1,2822√Re

Cabe comentar que si el intervalo en el eje X hubiera sido de 4x = 0,001 se habríaobtenido Cf = 1,263582√

Re, por lo que el error cometido se reduce por anamiento del mallado

como era de esperar.

En la gura 5.7 se observan los valores del coeciente de fricción teórico y del obtenidopor métodos numéricos para cada punto de la placa, así como el valor del coeciente defricción medio.

Este resultado se puede comparar con el teórico deducible de la teoría de Blasius.Retomando el valor teórico de Cf (x) = 0,664√

Ree integrando para toda la placa se halla el valor

medio de Cf :

Cf =

´ L0Cf (x)dx

L=

´ 10

0,664√Re√x′Ldx

L=

1,328√Re

46 ETSII-UPM


0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Longitud adimensionalizada de la placa

0

0.5

1

1.5

2

2.5

3

3.5

4

Cf(x)

Coeficientes de fricción teórico y numérico


Figura 5.7: Resultados obtenidos Cf en la placa plana

Por lo tanto, el valor obtenido es muy próximo al teórico, siendo el error cometido del3,4 %. Este error se debe a que en los primero puntos el perl de velocidades pasa de serigualmente a 1 para todos los nodos a un perl en el que el primer valor vale 0, y este cambiobrusco induce error en el primer 10 % de la placa plana. Además se comprueba que Cf esmayor que 0 y no nulo para todos los puntos, por lo que la capa límite no se desprende, comose había previsto.

5.3. Succión homogénea en la placa plana

El objetivo de este apartado no es más que presentar un método de control de la capalímite y de nuevo comparar con las soluciones que se presentan en la referencia [12], parajusticar la validez del método y vericar su correcto desarrollo.

La succión es un método que modica la capa límite de forma articial como se muestraen la gura 5.8. Consiste en la absorción de partículas uidas a través de una de las paredesdel cuerpo sobre el que se forma la capa límite. De esta manera la capa límite se hace de menorespesor, por lo que se retrasa su crecimiento y posible desprendimiento. En geometrías en lasque exista gradiente de presión adverso se retrasará el desprendimiento de la capa límite, sibien en el caso de la placa plana tan sólo se retrasará su crecimiento. Este método o similaresse utilizan habitualmente en las alas de los aviones, ya que se reduce la resistencia al avancepor disminución de la fuerza de presión Fp. Esto ocurre porque si se aplica succión se retrasala transición de la capa laminar a turbulenta, por lo que la resistencia al rozamiento es menor.



Figura 5.8: Succión en la capa límite

El caso que se va a estudiar es el de la succión homogénea, es decir, a lo largo de lageometría se extrae uido con una velocidad constante e igual a v0. Se va a aplicar la succiónal caso de la placa plana, así como al cilindro, como se verá en el siguiente capítulo.

5.3.1. Solución analítica

Para el caso de una placa plana con succión homogénea existe solución analítica [12].Las curvas que se van a obtener son las que aparecen en la gura 5.9.

Figura 5.9: Curvas para succión homogénea. Fuente [12]

La gura de la derecha muestra los perles de velocidades para distintos puntos de laplaca plana. De esta manera de representan los valores en función de −v0y

ν, donde v0 es la

velocidad de succión y ν = µρ. Por otro lado ξ es un parámetro tal que

√ξ = −vo

U∞

√U∞xν, y

representa distintos puntos x de la placa. Cabe comentar que en el programa correspondienteen Anexos la velocidad del succión v0 se ha denominado vs.

La gura de la izquierda se observan los perles de velocidades sin succión (perl II)y el perl de velocidades teórico que se obtiene si hay succión homogénea en la placa (perlI). Además y

δ= −v0y

ν.

48 ETSII-UPM


Se obtendrán estas curvas en Matlab, si bien antes se explica el signicado de cada unade ellas. El perl I tiene solución analítica, la cual se deduce a continuación.

Partiendo de las ecuaciones de Navier Stokes y sus condiciones de contorno:

u(y = 0) = 0

v(y = 0) = −v0u(y =∞) = U∞

∂u

∂x+∂v

∂y= 0

ρ(u∂u

∂x+ v

∂u

∂y) = ρUp

dUpdx

+ µ∂2u

∂y2

donde v0 es la velocidad de succión y Up la velocidad del ujo potencial, en este casoUp(y = 0) = 0. Particularizando en y = 0 se tiene:

v∂u

∂y|y=0 = µ

∂2u

∂y2|y=0

Integrando e imponiendo las condiciones de contorno se llega a que:

u(y) = U∞(1− eµv0y)

Este perl se llama perl asintótico de succión y según la bibliografía consultada [12] sealcanza en ξ = 4.

Con estos datos ya se puede decidir qué velocidad de succión tomar para resolver elproblema. Cabe comentar que se recomiendan velocidades del orden v0 ∼ −0,01U∞.

Antes de alcanzar el valor de v0 se van a adimensionalizar todos los parámetros expli-cados unas lineas más arriba para facilitar la comprensión del razonamiento. De esta manera,teniendo en cuenta:

x′ =x

L

y′ = y

√Re

L

u′ =u

U∞

v′ = v

√Re

U∞

Re =ρU∞L

µ



donde la variables con comillas son las adimensionalizadas. Para facilidad del lectorse retoma la notación anterior sin comillas para las variables adimensionalizadas. De estamanera se obtiene:

u(y) = 1− ev0y√ξ = −v0

√x (5.3)

y el eje horizontal de los grácos de la gura 5.9 es:

−v0y

Por lo tanto, si para ξ = 4 se alcanza el perl asintótico, y la longitud adimensionalizada dela placa es x = 1, se obtiene que la velocidad de succión adimensionalizada es v0 = −2 segúnla ecuación 5.3. Para vericar su valor se calcula cuánto vale:

v0U∞ =−2√Re

U∞ = −0,02U∞

Por lo que el orden de magnitud es correcto. Cabe preguntarse por qué la imposiciónde la geometría impone el valor de la velocidad de succión. Esto es porque el perl asintóticode succión se alcanza al nal de una placa de longitud innita. Como el problema que se haplanteado durante todo este capítulo se ha hecho para una longitud adimensionalizada de1, el nal de la placa de longitud real se alcanza en x = 1. Por lo tanto, se ha comprobadoque ξ = 4 para un punto de placa innito, es decir, sucientemente grande como para quese haya desarrollado la capa límite laminar, y que en este caso se ha considerado que es elpunto x = 1. Otra opción posible habría sido considerar el nal de la placa en, por ejemplo,x = 0,75, y se observaría que de x = 0,75 a x = 1 el perl de velocidades no varía, pues es elperl asintótico de succión.

A continuación se procede a la resolución del problema, para lo cual se han tenidoque cambiar las condiciones de contorno del problema de ujo potencial así como el de laresolución de la capa límite.

5.3.2. Flujo potencial

Para la resolución de este problema problema se ha tomado una geometría rectangularcomo la de la gura 5.2. Se han impuesto las condiciones de contorno tales que las funcionesde corriente recojan la información que provoca la succión, es decir, que las funciones decorriente muestran un caudal saliente por la placa plana y evidentemente en este caso lafunción de corriente no es constante sobre placa plana. Las dimensiones de la geometría síson relevantes en este caso, dado que la solución no es trivial. La altura del mallado ha de sersucientemente grande como para que la frontera superior se pueda suponer como una lineade corriente horizontal, y que no se ve afectada por la succión que experimenta el ujo en la

50 ETSII-UPM


frontera inferior. Se ha elegido de largo longitud unidad y de alto el triple. De esta maneralas funciones de corriente sobre cada frontera quedan:

ϕC1 = −wU∞ + v0 x

ϕC2 = Us(w − y)

ϕC3 = 0

ϕC4 = U∞(w − y)

donde Us es la velocidad de salida y las fronteras son las mismas que las de la gura 5.2.

Por lo tanto, para x = L, siendo L la longitud de la placa, que es la unidad, se tiene queϕC1 = −wU∞ + v0 L, que será igual a ϕC3. Se comprueba haciendo un balance de caudales:

wU∞ = v0 L+ wUs

Por lo que con este balance se verica la asociación de funciones de corriente a lasfronteras. Éstas pueden verse en la gura 5.10, donde se aprecia una desviación en las mismashacia la placa debido a la succión. Se recuerda que si no hubiera succión las funciones decorriente serían lineas horizontales paralelas a la placa.

Por lo tanto, una vez resuelto el campo de velocidades correspondiente al ujo poten-cial se tiene le vector de presiones a introducir en el programa de Matlab que resuelve lasecuaciones de la capa límite laminar.

5.3.3. Capa límite

En este caso hay que hacer un cambio en las condiciones de contorno respecto al caso dela solución de Blasius, que se recuerda que es la placa plana con ujo estacionario, continuoy paralelo a ésta. Las condiciones de contorno serán:

un0 = 0

vn0 = v0

unN = Unp = 1

u1(y) = u(x0, y)

es decir, en la formulación matricial del problema se tendrá que lN = v0√ReU∞

(4.8).

Se obtienen los perles de velocidades para este problema, que se comparan en la gura5.11. En este caso, se puede observar que el espesor de la capa límite es mucho menor, locual era de esperar, ya que como se comentó anteriormente, la capa límite se adhiere más ala placa.



Figura 5.10: Funciones de corriente en la placa plana para el caso de succión homogénea

En la gura 5.12 se han representado los perles de velocidades para los valores de√ξ = 0,1, 0,2, 0,4, 0,6, 0,8, 1 y 1,4, que como se puede comprobar, son iguales a los que se

muestran en la gura 5.9.

En la gura 5.13 se muestran tres perles. El primero de ellos es el perl obtenido sinsucción, el segundo el perl para ξ = 4, que según la bibliografía consultada debe coincidir conel perl teórico, y por último el perl asintótico de succión. Cabe comentar que para valoresde ξ superiores a 4 el perl obtenido sería el mismo. Sin embargo, como se ha explicado unaslíneas más arriba, se ha resuelto el problema hasta ξ = 4, que se corresponde con el perl dex = 1.

Como puede observarse el perl obtenido coincide con el teórico como se había esperado.Por otro lado, comparándolo con el perl que se obtiene de la solución de Blasius, el cual se

52 ETSII-UPM


0 0.5 1

u

0

1

2

3

4

5

6

7

Altu

ra

x=0.25

0 0.5 1

u

x=0.5

0 0.5 1

u

x=0.75

0 0.5 1

u

x=1


Figura 5.11: Comparación de los perles de velocidades.

obtuvo en la resolución del caso de la placa plana sin succión (gura 5.4), se observa que esmás lleno, pues la succión provoca una reducción en el espesor de la capa límite y el perlestá, consecuentemente, más curvado en la zona de la capa límite.

Finalmente, en este apartado se ha vericado la resolución del problema de succión enesta geometría.



0 0.5 1 1.5 2 2.5 3 3.5−yvsν

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

u Up

Perfil de velocidades en distintos puntos de la placa

√

ξ = 0.1√

ξ = 0.2√

ξ = 0.4√

ξ = 0.6√

ξ = 0.8√

ξ = 1√

ξ = 1.4

Figura 5.12: Distribución de velocidades en la placa placa con succión homogénea

0 0.5 1 1.5 2 2.5 3 3.5 4−yvsν

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

u Up

Comparación de los perfiles con y sin succión

Perfil sin succiónPerfil teórico con succiónPerfil obtenido con succión en el final de la placa

Figura 5.13: Distribución de velocidades en la placa placa con succión homogénea, sin succióny perl teórico

54 ETSII-UPM

Capítulo 6

Flujo alrededor de un cilindro

En este apartado se procede a analizar el campo uido alrededor de un cilindro segúnel método explicado: análisis del ujo potencial y posterior resolución de las ecuaciones de lacapa límite laminar con las condiciones de contorno de ujo potencial. Previamente se haráuna descripción cualitativa del cómo ha de ser el ujo alrededor del cilindro, pues es unageometría muy estudiada en la literatura.

6.1. Características del ujo

Sea una corriente estacionaria con velocidad U∞ = 1 incidiendo de izquierda a derechasobre un cilindro de longitud mucho mayor que diámetro. El punto situado más a la izquierdaes un punto de remanso, pues la velocidad del ujo es nula en ese punto. Tomando ese puntocomo origen de ángulos, para ángulos menores de 90, el gradiente de presiones es favorable,por lo que el ujo se acelera en ese tramo, alcanzándose la velocidad máxima del ujo en90. Para ángulo superiores a este último, el gradiente de presiones es adverso, por lo que elujo se frena existiendo la posibilidad de que se desprenda la capa límite formada, lo cualocurre si retrocede mucho el perl de velocidades. La gura 6.1 reeja este caso de manerailustrativa.

Para el caso capa límite turbulenta, que aparece con altos números de Reynolds, éstaresiste la separación de ujo mejor que un ujo laminar, ya que predominarían las fuerzasinerciales sobre las viscosas. El caso de capa límite turbulenta se sale del temario abarcadoen este proyecto, si bien un ejemplo muy ilustrativo es la inuencia del estado supercialdel cuerpo sobre el carácter de la capa límite. Así, si la rugosidad supercial es elevada, sefavorece la transición de capa límite laminar a turbulenta, por lo que se separa la capa límitemás tarde. Este fenómeno es el que ocurre es las pelotas de golf, como se ilustra en la gura6.2.

55

CAPÍTULO 6. FLUJO ALREDEDOR DE UN CILINDRO

ϕ

∂p∂x > 0 ∂p

∂x < 0

Figura 6.1: Gradiente de presión alrededor de un cilindro

Por lo tanto, en este apartado se va a estudiar un cilindro, teniendo en cuenta que lasecuaciones de la capa límite son para régimen laminar, y que el ujo potencial correspondea números de Reynolds innitos, se busca un número de Reynolds de compromiso: 104.

Según lo visto hasta ahora, se espera un ángulo de desprendimiento de más de 90 alser el perl de presiones obtenido el correspondiente un campo de ujo potencial.

Para el problema de ujo potencial exterior a la capa límite existe solución analítica.Ésta fue obtenida por Schlichting [12] y predice el punto de desprendimiento de la capalímite en 108,8. A este resultado se llega resolviendo la ecuación que da la solución analítica,al igual que se hizo en el caso de la placa plana. Para el caso del cilindro, la resoluciónes de complejidad mucho mayor, por lo que la comparación de resultados será meramentecualitativa como se verá más adelante. Otro dato que se ha tenido en cuenta para resolver esteproblema, también tomado de la resolución de Schlichting, es la distribución de velocidadesalrededor de un cilindro de radio R para ujo irrotacional, que es la siguiente:

Up(x) = 2U∞senx

R(6.1)

Como el método desarrollado por Schlichting es el mismo que el que se pretende resolveren este proyecto, será un requisito llegar a una solución con punto de desprendimiento similara este.

56 ETSII-UPM


Figura 6.2: Inuencia de la rugosidad en el carácter de la capa límite

Por otro lado, como el problema es simétrico se ha resuelto la mitad superior delcilindro. Esta hipótesis se puede realizar porque se estudia el punto de desprendimiento de lacapa límite, no lo que ocurre una vez desprendida, dado que en estado estacionario dejaría deser simétrico el ujo y además se recuerda que este método deja de ser válido para el estudiodel ujo a partir del punto de desprendimiento.

En cuanto a la geometría, se ha tomado un cilindro de radio R, que será el parámetroutilizado para adimensionalizar las variables, de tal manera que

Re =ρU∞R

µ

Cabe comentar que la solución no depende del radio elegido como se comenta más adelante,tan solo depende del número de componentes del mallado espacial. De hecho se han probadootros radios y las leves diferencias observadas se deben a la diferencia de número de compo-nentes en vector que recorre la supercie de la semicircunferencia y la consecuente variacióndel resultado por variación del paso en el mallado.

6.2. Flujo potencial

La geometría tomada para resolver el ujo potencial ha sido la que se muestra en laparte izquierda de la gura 6.3, de tal manera que w = 4R. Se ha tomando radio R = 3,



que es el radio adimensionalizado con R, por lo que el radio real medirá un tercio del queaparece en el mallado de FreeFem++. Otra posible interpretación de esta dimensión es tomarla geometría de FreeFem++ como si tuviera dimensiones, es decir, el radio fuera R = 3m.Esto no es un problema, pues al pasar el vector de presiones al programa de Matlab seadimensionalizará como se explica más adelante.

Por otro lado, en la gura 6.3 se puede ver la geometría en la que se va a resolver elujo potencial y en la que se nombran las seis fronteras del dominio, donde se han impuestolas siguientes condiciones de contorno:

ϕC1 = ϕC3 = U∞y

ϕC2 = 3wU∞

ϕC4 = ϕC5 = ϕC6 = 0

La altura del dominio es elevada para que las funciones de corriente en el extremosuperior sean constantes como se ha impuesto. En esta zona el ujo no se ve afectado por laexistencia de un cilindro si está sucientemente alejado según la teoría del ujo potencial.

En cuanto al mallado, se ha aumentado el número de elementos en la zona del cilindroya que es la que interesa de cara a la resolución del problema de la capa límite.

La resolución del problema ha dado las funciones de corriente de la parte derecha de lagura 6.3 donde se observa que la función de corriente en el extremo superior es prácticamenteconstante como se pretendía.

El gradiente de presiones obtenido es el que se muestra en la parte derecha de la gura6.5. Se puede observar que es negativo hasta π

2, como era de esperar según lo comentado

al inicio de este capítulo. Por el contrario, a partir de π2el gradiente de presiones se hace

positivo, por lo que se tiene un gradiente adverso de presiones que será el que provocará eldesprendimiento de la capa límite en un modelo tan simplicado como este. Cabe comentarque este gradiente de presiones no se ha medido sobre la supercie de la circunferencia, sinoh = αR por encima, donde α se ha ido variando hasta obtener la solución óptima. Hay dosrazones para haber realizado esto. La primera de ellas es que FreeFem++ construye el cilindrode tal manera que no es una circunferencia exacta, sino una serie de nodos muy juntos unidospor rectas, y que todos ellos forman una semicircunferencia. Es por ello que la lectura de lapresión ahí conlleva mucho más error que en una zona donde hay más nodos alrededor, comopueda ser alejándose un poco más de la circunferencia. La segunda razón es coger una alturarealista que se corresponda con la zona exterior de la capa límite.

La solución óptima o correcta es aquella que hace el coeciente de presión lo más pare-cido posible al que resulta de la teoría de uidos no viscosos, lo cual se explica a continuación.

El caso que se está estudiando es un cilindro, por lo tanto el cuerpo es totalmentesimétrico respecto al plano horizontal. Consecuentemente el sumatorio de fuerzas en sentido

58 ETSII-UPM


Figura 6.3: Mallado sobre el cilindro y funciones de corriente

vertical es nulo al anularse por simetría geométrica. De esta manera la única fuerza que quedapara estudiar es la fuerza de arrastre D, que deriva de las fuerzas ejercidas por la presióny por los esfuerzos cortantes en el sentido horizontal. Adimensionalizando la fuerza ejercidapor la presión se obtiene el coeciente de presión en cada punto del cilindro, que se denecomo:

Cp(x) =(p− p∞)

12ρU2∞

Con las variables adimensionales que se han tomado a lo largo del proyecto queda:

Cp(x) =ρU2∞p(x)

12ρU2∞

= 2p(x)

Con esta fórmula el coeciente de presión se puede calcular en cada punto si las presionesson conocidas, las cuales se obtienen de la resolución del campo potencial.

Por otro lado, el coeciente de presión teórico tiene la siguiente forma [14]:

Cp(x) = 1− 4sen2(ϕ)



0 0.5 1 1.5 2 2.5 3 3.5α

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

Cp(x)

Coeficiente de presión Cp en función de α

h = 0.01R

h = 0.02R

h = 0.03R

h = 0.05R

h = 0.1R

Cp teórico

Figura 6.4: Valores de Cp en función de h

Por lo tanto, probando para varios valores de h = αR para una radio de R = 3 se obtiene lagráca que se puede ver en la gura 6.4.

En vista de los valores obtenidos, se elige α = 0,025, ya que es el valor que se encuentraentre las curvas de Cp para h = 0,02R y h = 0,03R. También se pueden ver los gradientesde presiones para los distintos valores de h en la gura 6.5. Finalmente, para h = 0,025Rel gradiente de presiones obtenido es el que se muestra en la parte derecha de la gura 6.5.Por el razonamiento seguido, el gradiente de presiones será el mismo que el teórico, lo cualverica los cálculos.

El programa de FreeFem++ almacena un valor de la presión cada 0,01 de semicircun-ferencia recorrida. De esta manera, el vector de presiones tendrá πR

0,01elementos.

6.3. Capa límite

El vector x que recorre la semicircunferencia se ha adimensionalizado con R, por lo queirá de 0 a π. El vector y′′ =

√ReRy va de 0 a 3,5, por ser este el rango óptimo que abarca la

evolución del perl de velocidades en la capa límite. Si fuera hasta más valores superiores a3,5 las ecuaciones resolverían parte del campo de velocidades que está fuera de la capa límite,de tal manera que para los nodos superiores las velocidades serían constantes e iguales a ladel ujo potencial, pues no son válidas las ecuaciones de la capa límite en esta región. Es

60 ETSII-UPM


0 0.5 1 1.5 2 2.5 3

Longitud de la circunferencia

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

2.5Gradiente de presiones

h = 0.01R

h = 0.02R

h = 0.03R

h = 0.05R

h = 0.1R

0 0.5 1 1.5 2 2.5 3 3.5

Longitud de la circunferencia adimensionalizada

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

2.5Gradiente de presiones

Figura 6.5: Gradiente de presiones a lo largo del cilindro

decir, si se resuelve más dominio que el que abarca el espesor de la capa límite se impondráuna velocidad igual a la del ujo potencial hasta el innito, lo cual, evidentemente, no esreal. Si el valor fuera inferior la resolución no recogería toda la capa límite y sería errónea. Denuevo, este paso en el mallado se ha determinado probando distintos y evaluando la validezde los perles de la capa límite. Al igual que para otras geometrías, para esta se ha elegido4y = 0,01.

En número de elementos del vector x es el mismo que el del vector de presiones queda FreeFem++, y a partir de eso se obtiene el valor de 4x. De esta manera, el númerode componentes de ambos vectores es Nx. Esta es la explicación de que no importa si lageometría de FreeFem++ esté con o sin dimensiones, pues lo único que importa es el númerode componentes del vector de presiones:

Nx =πR

0,01

Nx =π

4x

Por lo tanto 4x = πNx, que es el paso del mallado de Matlab en sentido horizontal

(sección 4.5.1).

Cabe comentar que se han probado otras longitudes del vector y y otros valores de 4xy 4y, siendo nalmente los elegidos los que recogen correctamente la evolución del perl alo largo de la semicircunferencia sin necesidad de mayor exactitud.

En este caso la condición de contorno a imponer en los nodos UnN para cada n es la que

se ha obtenido de la resolución de Schlichting, por lo que la ecuación 6.1 teniendo en cuentalas variables adimensionalizadas queda como:

UnN = Un

p = 2U∞ sen x(n)



donde se recuerda que N es el nodo superior de cada columna.

Resolviendo el problema se han obtenido los perles de velocidades en la capa límite,y se representan tres de ellos en la gura 6.6. Los puntos representados son ϕ = π

4, ϕ = π

2y

el punto de desprendimiento obtenido.

0 0.5 1 1.5

u

0

0.5

1

1.5

2

2.5

3

3.5

Altu

ra

ϕ =π

4

0 1 2

u

ϕ =π

2

0 1 2

u

ϕ = 105.48

Figura 6.6: Perles de velocidades en el cilindro

-0.03 -0.02 -0.01 0 0.01 0.02 0.03

u

ϕ = 105.48

Figura 6.7: Ampliación del último perl de velocidades en el punto de desprendimiento.

La longitud de la semicircunferencia adimensionalizada es π. El programa resuelve hasta1,84, punto que se corresponde con el de desprendimiento de la capa límite. Como se puedever en el código del programa, se detiene el cálculo cuando el perl de velocidades se da lavuelta, es decir, se hace negativo en la parte más inferior. Si se amplía la parte inferior del

62 ETSII-UPM


perl correspondiente al punto de desprendimiento se puede apreciar el perl dado la vueltacomo se muestra en la gura 6.7. Ese punto de desprendimiento se corresponde con un ángulode desprendimiento

ϕ = 105,4777

Conocidos los perles de velocidades, se procede a compararlos con los resultantes de lasolución analítica de la capa límite. En la gura 6.8 se presentan los perles de velocidades devarios ángulos. En el eje horizontal se representa la velocidad horizontal u adimensionalizadacon Un

p .

0 0.2 0.4 0.6 0.8 1uU

0

0.5

1

1.5

2

2.5

3

y R

√

URν

ϕ = 20

ϕ = 40

ϕ = 60

ϕ = 80

ϕ = 90

ϕ = 100

ϕ = 105.48

Figura 6.8: Comparación de las curvas obtenidas. Fuente [12]

Como se puede observar, las curvas son prácticamente iguales. Las diferencias quehay se deben a la inexactitud del método, motivada por el paso del vector de presiones deFreeFem++ a Matlab y por el uso de un número de Reynolds de teoría laminar para unvector de presiones de ujo potencial. Estas inexactitudes también se ven reejadas en queel ángulo de desprendimiento no sea exactamente 108,8, y el error cometido es del 3,05 %.

Finalmente se muestra el coeciente de presiones teórico y el obtenido, que como ya sehabía adelantado, coincide con el teórico, pues a partir de esa condición se ha determinadoh = αR. Todo ello se puede ver en la gura 6.9.



0 0.5 1 1.5 2 2.5 3 3.5ϕ

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

Cp(x)

Coeficiente de presión en función de la posición

Coeficiente de presión obtenidoCoeficiente de presión teórico

Figura 6.9: Coeciente de presiones en función de la posición

Cabe comentar que el coeciente de arrastre CD para el cilindro en función del númerode Reynolds está tabulado. Sin embargo, el cálculo del mismo en este caso no va a resultarcorrecto, puesto que el programa realizado resuelve las ecuaciones de la capa límite hastael punto de desprendimiento y a partir de éste, deja de ser válida esta teoría. El coecientede arrastre que otorga la solución de las ecuaciones de Navier Stokes y que está tabuladoen la mayoría de la bibliografía consultada está calculado para el cilindro completo, inclusodespués del desprendimiento de la capa límite, por lo que no tendría sentido compararlos.

Por lo tanto, en vista de todos los resultados obtenidos y contrastados, de nuevo seconrma la validez de este método.

6.4. Succión homogénea en el cilindro

Como se ha visto anteriormente, la succión es una manera de retrasar el desprendimien-to de la capa límite. En este caso no hay solución analítica pero se explicará el procedimientode cálculo de los perles de velocidades en el caso de succión homogénea a lo largo del cilindro.Finalmente se presentará la mejora conseguida de forma cualitativa.

Se ha decidido resolver con una velocidad de succión similar a la utilizada para el casode la placa plana, es decir, del orden de 0,01U∞.

64 ETSII-UPM



Para la resolución de este problema se ha tomado un mallado como el de la gura6.3. Se han impuesto las condiciones de contorno tales que las funciones de corriente recojanla información que provoca la succión, es decir, que la funciones de corriente muestran uncaudal saliente por la semicircunferencia. De esta manera:

ϕC1 = Us(y − 3w)

ϕC2 = 0

ϕC3 = U∞(y − 3w)

ϕC4 = −U∞3w

ϕC5 = −U∞3w + v0Rθ

donde Us es la velocidad de salida y v0 la velocidad de succión. Además θ es el ángulomedido desde el punto de remanso, que hasta ahora se ha nombrado como ϕ pero en esteapartado se cambia la notación para evitar confusiones con las lineas de corriente. A diferenciade la resolución en el programa sin succión, la función de corriente nula se ha impuesto sobrela frontera C2, lo cual se ha hecho por comodidad. Además se ha elegido v0 = 0,02U∞.

Por lo tanto, para x = πR, siendo R el radio del cilindro, se tiene que ϕC3 = −3wU∞+v0Rθ, que será igual a ϕC3. Se comprueba haciendo un balance de caudales:

3wU∞ = v0 πR + 3wUs

Las líneas de corriente pueden verse en la gura 6.10, y se aprecia una desviación en lasmismas hacia la placa debido a la succión. Si no hubiera succión las funciones de corrienteserían las mostradas en la gura 6.3.

6.4.2. Capa límite

En este caso hay que hacer un cambio en las condiciones de contorno respecto al casosin succión. Las condiciones de contorno serán:

un1 = 0

vn1 = v0

unN = Unp = 1

u1(y) = u(x0, y)

es decir, en la formulación matricial del problema se tendrá que lN = v0√ReU∞

.



Figura 6.10: Funciones de corriente con succión homogénea en el cilindro

Los perles de velocidades obtenidos se muestran en la gura 6.11, así como lo que seobtienen si no hay succión. Se puede observar una mejora en el espesor de la capa límite, esdecir, ésta se hace más estrecha cuando hay succión. Sin embargo, la mejora más notable esel retraso del punto de desprendimiento. La gura 6.11 resulta de resolver el problema conv0 = 0,02. Si bien en la tabla 6.1 se muestran los distintos ángulos de desprendimiento enfunción de la velocidad de succión.

Velocidad de0,00 0,001 0,005 0,01 0,015 0,02 0,025 0,03

succión v0U∞

Ángulo de105,478 106,815 112,548 120,764 129,554 138,917 148,471 158,026

desprendimiento ϕ

Cuadro 6.1: Variación del ángulo de desprendimiento según la velocidad de succión

Las conclusiones de este apartado de succión son claras. Por un lado la succión retrasael punto de desprendimiento por con respecto al que se obtiene sin succión. Además, a mayorvelocidad de succión, más retrasado está el punto de desprendimiento, lo cual se aprecia enla tabla.

Como conclusión del capítulo, se verica de nuevo la validez del método, pasando a laresolución de la aplicación del mismo en el próximo capítulo.

66 ETSII-UPM


0 0.5 1 1.5

u

0

0.5

1

1.5

2

2.5

3

3.5

Altu

ra

ϕ =π

4


Figura 6.11: Comparación de los perles de velocidades



68 ETSII-UPM

Capítulo 7

Aplicación: ujo en arterias con estenosis

Como aplicación del modelo matemático desarrollado se va a simular el ujo sanguíneoen una geometría especíca para su posterior análisis. La aplicación estudiada es la estenosisarterial, una enfermedad que consiste en una variación de la geometría de los vasos sanguíneosque afecta al ujo y cambia sus condiciones locales del mismo. La estenosis es un términoque denota el estrechamiento de un oricio o conducto corporal. Algunos de los tipos máscomunes son la estenosis de espinal, intestinal o aórtica.

La estenosis arterial, concretamente, es un estrechamiento del vaso sanguíneo que pue-de provocar el desprendimiento de la capa límite, produciéndose recirculaciones que llevana la deposición de lípidos. Este estrechamiento puede llegar a ocluir totalmente la luz delvaso. La causa más frecuente de la estenosis la aterosclerosis, que es una enfermedad ina-matoria crónica de las arterias que consiste en el engrosado progresivo como consecuenciade la acumulación de diversas sustancias, ya sean grasas, calcio o células de la sangre. Elengrosamiento de la pared arterial provoca una dicultad al paso de la circulación sanguínea,pudiendo llegar a obstruir completamente la arteria. Estos engrosamientos se conocen comoplacas de ateroma. Generalmente, estas patologías aparecen en bifurcaciones de la aorta (bi-furcación carotídea, en el cuello) o arterias renales, así como en arterias más pequeñas. Sedetectan a través de angiografías, y en la gura 7.1 se puede apreciar la geometría real deesta enfermedad.

En cuanto al tratamiento médico de esta enfermedad, es necesario un control rigurososobre los factores de riesgo cardiovascular implicados en el crecimiento de la placa de ateroma,así como la toma de fármacos. Además existe otra forma de tratamiento más agresivo queconsiste en la apertura forzada de la luz de la arteria mediante cirugía directa sobre el vasoo a través de un catéter introducido por la ingle para colocar un stent.

La ampliamente conocida trombosis está relacionada con fenómenos uidodinámicoscomo puedan ser altos valores del esfuerzo cortante en la pared del vaso sanguíneo, recircula-ciones, separación de la capa límite o puntos de remanso. Este tipo de anomalías geométricaspueden llevar a la migración de glóbulos rojos alejándose de la pared, por lo que la concen-tración de plaquetas aumenta en la pared.

69

CAPÍTULO 7. APLICACIÓN: FLUJO EN ARTERIAS CON ESTENOSIS

Figura 7.1: Angiografía de arteria carótida. Estrechamiento del 80% [6]

Se han realizado numerosos experimentos cuyo objetivo es ayudar a la medicina apredecir el comportamiento de la sangre en estos casos. Estos experimentos se llevan a cabohabitualmente con muestras de sangre de animales como la oveja, el mono o el perro, cuyosparámetros son parecidos a los de la sangre humana.

A partir del punto de desprendimiento de la capa límite el método utilizado no puedepredecir el comportamiento del ujo, si bien cabe comentar que en muchos casos la capalímite se desprende y, tras recorrer cierta distancia, ésta se reengancha. Por lo general ladistancia entre los puntos de desprendimiento y reenganche es proporcional al número deReynolds, lo cual se comprueba de forma experimental.

También se observa experimentalmente que el ujo constante y pulsátil no presentalesiones debidas a un alto esfuerzo cortante en la pared del vaso [9]. Sin embargo, el ujoinestable o no estacionario presenta regiones de recirculación que están relacionadas con laszonas donde se producen las lesiones. Como ya se ha comentado, las arterias que sufrenestenosis producen lesiones por acumulación de depósitos en la pared arterial. Por lo tantoes de interés cientíco el estudio del ujo no estacionario, si bien en este proyecto se hará unestudio del ujo estacionario al igual que las geometrías presentadas hasta ahora.

Se han estudiado distintas geometrías de la enfermedad y para cada una de ellas,distintas condiciones del ujo. Según la Organización Mundial de la Salud existen tres gradosde severidad en la enfermedad en función de la reducción de la luz de la arteria:

Sin estenosis: si apenas reducción de la luz arterial.

Estenosis leve: reducción del lumen de arteria de menos del 54%.

Estenosis moderada: si todavía la luz de la arteria en el estrechamiento es mayor de lamitad del diámetro del lumen original, es decir, reducciones de hasta el 75%.

70 ETSII-UPM


Estenosis severa: si la reducción de la luz arterial es mayor del 75%, o lo que es lomismo, el diámetro en el estrechamiento es menor que la mitad del original.

Aunque hubo ciertas discrepancias hace años, hoy en día, a la vista de numerosos ensayos,se considera que las zonas más proclives a la estenosis son las sometidas a bajo esfuerzocortante, como en las recirculaciones, donde se depositan las placas de ateroma. El depósitoformado provoca el engrosamiento de la pared o estenosis con la consiguiente aceleracióndel ujo y aumento del esfuerzo cortante. Éste último es el responsable de la activaciónplaquetaria y trombosis. Por lo tanto se estudiará el valor del esfuerzo cortante en las paredesde vasos sanguíneos con estenosis para posteriormente comparar los resultados obtenidos conlos experimentales. Es decir, el alto esfuerzo cortante se toma como indicador de la gravedadde la estenosis y consecuentemente se estudia la enfermedad a través de este parámetro, todoello en vasos sanguíneos que ya tienen estenosis.

Cabe comentar que el máximo valor del esfuerzo cortante se da antes del estrecha-miento máximo del vaso y, como es predecible, es mayor cuanto mayor es el estrechamiento.Se comprueba experimentalmente que si este parámetro supera el valor de 380 dina cm−2

puede aparecer disfunción endotelial [9]. Esta es una de las primeras manifestaciones de laenfermedad vascular y la arteriosclerosis. El endotelio, una monocapa de células que recubrela pared luminal de los vasos sanguíneos, regula la interacción de las células y las proteínascirculantes con las células residentes en la pared vascular, ejerciendo un papel central comosensor y transmisor de señales.

De esta manera, queda claro que un alto esfuerzo cortante en la pared puede provocardisfunción endotelial, con el consiguiente inicio de la arteriosclerosis. La estenosis puede serconsecuencia de ésta, y llevar a la obstrucción de la arteria por formación de un coágulo,entre otras cosas.

Por lo tanto, el estudio del esfuerzo cortante para distintas geometrías y distintos ujospermite analizar la etapa en la que se encuentra la enfermedad. Cabe comentar que enel artículo de Huang [9] se arma que un alto esfuerzo cortante puede no ser causa de laarteriosclerosis de forma cualitativa, ya que en ciertas referencias se encontró que no habíalesión la parte de la contracción y sí en la parte distal y proximal. Además en esos animales nose encontró disminución del ujo para contracciones de más del 75 %, por lo que esto podríanser resultados conservadores y, de forma cualitativa, podría decirse que un alto esfuerzocortante no es causa de arteriosclerosis.

7.1. Geometría

Para realizar el estudio por métodos numéricos se ha empleado la geometría de la gura7.2, donde L es la longitud de la estenosis y H la altura de la misma. La reducción será:

A = 1− R22

R21



Esta geometría cónica es una aproximación de la geometría real de la estenosis, quepor supuesto podrá ser mejorada, lo cual que verá más adelante.

R1

L

R2

H

Figura 7.2: Geometría de la estenosis

Los distintos casos estudiados en la bibliografía [9, 10] son los que aparecen en elcuadro 7.1. Estos casos se corresponden con los distintos grados de severidad de la estenosis:leve(M1), moderada (M2) y severa (M3, M4 y M5).

Parámetros M1 M2 M3 M4 M5HR1

14

13

12

12

12

LR1

4 4 4 8 2

A( %) 44 56 75 75 75

Cuadro 7.1: Parámetros geométricos para el modelo

7.2. Características del ujo

El ujo sanguíneo se considera laminar hasta Re = 2000, y turbulento a partir deRe = 3000. Entre estos dos valores hay una transición entre ujo laminar y turbulento.

El modelo matemático desarrollado se aplicará para el estudio del ujo sanguíneo enarterias con estenosis para números de Reynolds entre 100 y 1000, por lo tanto ujo laminar.

Se asume que el ujo es estacionario y no pulsátil, así como una geometría simplicadade la estenosis, pues en la realidad es irregular. Con el mismo n, se simplica el cálculosuponiendo que el ujo está completamente desarrollado al inicio del estrechamiento, y quese trata de un ujo regido por la ecuación del ujo de Poiseuille. Como se explica más adelante,

72 ETSII-UPM


esta geometría se resuelve por dos métodos distintos. El ujo de Poiseuille a la entrada setoma para uno de ellos, si bien para el otro será un perl recto de velocidad uniforme. Eluido se considera newtoniano y las paredes rígidas, es decir, se desprecia la elasticidad delos vasos sanguíneos. El uido es incompresible y su viscosidad es constante.

Los parámetros necesarios para conocer las características del ujo y de la geometría sonlos siguientes: ρ = 1060 kg m−3, µ = 3,71 10−3N sm−2 y R1 = 3,5 10−3mm [9, 10, 4]. Estosdatos se miden experimentalmente, para luego poder comparar los resultados obtenidos conlos experimentales. El radio podría corresponder al de una arteria como pueda ser la carótidainterna, situada en el cuello [5].

En cuanto a la viscosidad, antes de la estenosis es signicativa ya que el perl de velo-cidades corresponderá al ujo de Poiseuille en uno de los métodos, y está presente en toda laanchura de la arteria. En cambio, en el estrechamiento se ha supuesto un perl de velocidadesdistinto como se sugiere en el artículo de Reese [10]. De esta manera, en el estrechamientola viscosidad solo será importante en una pequeña región, donde está la capa límite. En lagura 7.3 pueden verse ambos perles.

Figura 7.3: Perles de velocidades en la estenosis

La velocidad a la entrada presenta un perl de Poiseuille, cuya velocidad máxima esU1, mientras que el estrechamiento será U2, y el perl de velocidades no es el del ujo dePoiseuille.

7.3. Datos experimentales

Los datos que se presentan en el cuadro 7.2 corresponden al esfuerzo cortante máximomedido experimentalmente en animales, para el caso de ujo de Poiseuille completamentedesarrollado. Se compararán con éstos los resultados obtenidos. Estos datos se han obtenidoa través de modelos experimentales que simulan el ujo sanguíneo a través de la estenosis[10].



Re M1 M2 M3 M4 M5

100 8,3 13,2 36,7 28,6 49,9500 67,6 110,6 316,9 240,3 423,21000 168,3 279,2 788,1 607,7 1022,7

Cuadro 7.2: Esfuerzo cortante máximo experimental (dina cm−2)

7.4. Metodología

Se va a resolver el problema de una manera ligeramente distinta a la que se ha realizadohasta ahora. En resumen, no se va a aplicar teoría del ujo potencial a la parte exterior dela capa límite, pues para números de Reynolds bajos, como es sabido, no se debe aplicar.En este caso, por lo tanto, las condiciones de contorno de las ecuaciones de la capa límitelaminar están relacionadas con el hecho de que el ujo sanguíneo está regido por la ecuaciónde Poiseuille en la zona exterior a la capa límite.

Teniendo en cuenta que en un vaso sanguíneo se tiene un ujo regido por la ecuaciónde Poiseuille salvo la zona de la capa límite, se calculará el esfuerzo cortante máximo y secomparará con el experimental.

Además, para relacionarlo con el ujo potencial, se resolverá está misma geometríapara el caso de números de Reynolds altos (Re = 2000), por lo que la aproximación ala teoría del ujo potencial será correcta, como se demostrará. Para este caso se utilizará lametodología desarrollada en los capítulos anteriores: imposición de condiciones de contorno deujo potencial a las ecuaciones de la capa límite laminar. La comparación entre los resultadospara ujos de Re = 2000 con condiciones de contorno de teoría potencial y de ujo dePoiseuille se llevará a cabo a través del estudio de la posición de los puntos de desprendimientode la capa límite. Se recuerda que la teoría potencial es aplicable a números de Reynolds muyaltos, tendiendo a innito, por lo que el uido se aproxima a uido ideal (sin viscosidad) fuerade la capa límite. El acoplamiento de la solución de ujo potencial a las ecuaciones de la capalímite laminar requiere que se trabaje con números de Reynolds lo más altos posibles perolaminares, como se ha hecho en los capítulos anteriores. Por otro lado, la teoría de Poiseuillese aplica a uidos en los que la viscosidad no es despreciable, como es el caso de la sangre. Porlo tanto, el desprendimiento de la capa límite se retrasará para el caso del ujo de Poiseuilleen comparación con los que se obtienen de la resolución de las ecuaciones de la capa límitelaminar, pues los esfuerzos viscosos son los responsables de la adherencia de la capa límiteal cuerpo. Consecuentemente, se puede adelantar que los puntos de desprendimiento seránanteriores para el método de resolución ampliado a Re = 2000 con ujo potencial.

74 ETSII-UPM


7.5. Resolución del ujo sanguíneo con condiciones de

contorno de ujo de Poiseuille

La resolución de este apartado se explicará en las siguientes líneas, y se compararánlos resultados obtenidos con los experimentales y con obtenidos por otros autores por reso-lución de las ecuaciones de Navier Stokes u otros métodos y que se muestran en los artículosconsultados. Todo ello se ha resuelto en Matlab.

Para implantar el modelo matemático se resolverán las ecuaciones de la capa lími-te imponiendo las condiciones de contorno siguientes: perl inicial de Poiseuille conoci-do y velocidad en la garganta también. Los datos del problema son ρ = 1060 kg m−3,µ = 3,71 10−3N sm−2, R1 = 3,5 10−3mm, Re, A y L, que dependerán del caso a estudiar.Las variables que se utilizarán para adimansionalizar serán R1 y U1 = µRe

ρR1.. Por otro lado

R2 = R1

√1− A.

El perl inicial será de Poiseuille como ya se ha comentado:

up(r) =U1

R2(r2 −R2)

Este perl es parabólico, y r es la distancia radial medida desde el eje de la tubería deradio R, por lo que se hace nula la velocidad en las paredes del vaso y máxima en el eje delmismo. Teniendo en cuenta las adimensionalizaciones del problema y′′ =

√ReR1y, u′ = u

U1y que

r2 = (R− y)2, el perl adimensionalizado queda:

u′p =y′′√Re

(2− y′′√Re

) (7.1)

Por lo que u′p = 1 si y = R1. Por otro lado la velocidad media, según el ujo dePoiseuille, es Um = U1

2. Por lo tanto, por conservación del caudal a lo largo del conducto se

tiene:

Q = πR21

U1

2= πR2

2U2 (7.2)

U2 =U1R

21

2R22

7.5.1. Corrección de la velocidad con el espesor de desplazamiento

Como es de esperar el ujo se acelera en la garganta por la disminución de la luz arterial.Sin embargo, la ecuación 7.2 debe ser corregida restando es espesor de desplazamiento de lacapa límite δ∗ a R2. A continuación se explica la razón.



En la gura 7.4 se muestra en color azul el perl real que se obtendrá en la garganta,desde la pared del vaso hasta el eje, por lo que la altura es R2, y en color rojo el que se estásuponiendo para calcular la velocidad U2. Como el espesor de desplazamiento no es conocido,se ha decidido restar una cantidad dR2 y se demuestra a continuación que es el espesor dedesplazamiento, que dene como:

δ∗ =

ˆ ∞0

(1− u

U) dy '

ˆ δ

0

(1− u

U) dy (7.3)

donde δ es el espesor de la capa límite, u la velocidad horizontal en la capa límite, y Ula velocidad exterior a ésta, es decir, en este caso U2.

Por lo tanto, para que sea válido, estos dos caudales han de ser iguales, es decir, en dosdimensiones se tiene: ˆ R2

dR2

U2dy =

ˆ R2

0

u dy

dR2 = R2−ˆ R2

0

u

U2

dy =

ˆ R2

0

(1− u

U2

) dy =

ˆ δ

0

(1− u

U2

) dy+

ˆ R2

δ

(1−U2

U2

) dy =

ˆ δ

0

(1− u

U2

) dy

Como se puede comprobar, esta última es la denición de espesor de desplazamiento(7.3), por lo que

dR2 = δ∗

R2

dR2

Figura 7.4: Espesor de desplazamiento

De esta manera se tiene la siguiente equivalencia:

76 ETSII-UPM


U2 =U1R

21

2(R2 − δ)2=

U1R21

2(R2 − dR2)2

El factor de corrección d se ha elegido inicialmente de forma arbitraria midiendo sobrelos perles de velocidad experimentales que se presentan en el artículo Huang [9]. Está claroque a mayor número de Reynolds el espesor de la capa límite es menor, por lo que d irávariando para cada número de Reynolds. A partir de ahí se ha iniciado un proceso iterativoteniendo en cuenta los datos experimentales, y nalmente se han obtenido valores de d paracada número de Reynolds.

Cabe comentar que para el caso de una placa plana el espesor de la capa límite varíade la siguiente manera, como ya se vio anteriormente (3.14).

δ

x' 5

Re1/2x

(7.4)

para régimen laminar yδ

x' 0,16

Re1/7x

para régimen turbulento.

Los valores obtenidos para los tres números de Reynolds para los que existen datosexperimentales son los que se muestra en el cuadro 7.3.

Re d

100 0,18500 0,101000 0,05

Cuadro 7.3: Parámetro de corrección d para los distintos Re

Con estos valores se ha determinado una ley de la forma

δ

R1

' A

Re1/Bx

donde A y B son constantes a determinar.

Lo ideal es obtener una ley similar a la de la placa plana (7.4), especialmente que latendencia sea la misma, es decir, B ' 2.

dα1

Re1/Bx

(7.5)

Para saber qué aspecto tiene esta ley deducida experimentalmente tras un proceso deiteración y comparación con los valores solución experimentales, se presenta en la gura 7.5



la curva de la ecuación 7.4 y la curva de una ecuación proporcional a 7.5. Estas curvas norepresentan el mismo parámetro, pues en rojo está representado el espesor de la capa límiteadimensionalizada con la longitud de la placa L para el caso de ujo laminar alrededor deuna placa plana y en azul está representado d, que es el espesor de desplazamiento de lacapa límite adimensionalizado con el radio del vaso R2. Lo relevante de esta gráca es quelas rectas son paralelas, es decir, B ' 2, de hecho, se ha calculado la ley 7.5 ajustando pormínimos cuadrados y se ha obtenido B = 1,9133.

101 102 103 104

Re

10-2

10-1

100

101

d

Ley δ=f(Re) para flujo laminar

Flujo sanguineo

Placa plana

Figura 7.5: Variación de d con respecto a la variación de δ en el caso de una placa plana

Otra de las consecuencias que supone que B sea aproximadamente 2 es la validez delcambio de variable que se hizo inicialmente y′′ =

√ReR1y, es decir, las variables de espacio ver-

tical adimensionalizadas son proporcionales a la raíz del número de Reynolds por la variablesin adimensionalizar: y′′α

√Rey. Por lo tanto, lo que se ha obtenido es que d

√Re = cte, es

decir, δ∗√Re = cte.

Cabe comentar que esta ley es válida para régimen laminar, al igual que la que imperaen la placa plana (7.4), lo cual ocurría cuando Re ' 2000, como ya se ha comentado. Esteparámetro se utilizará más adelante para comparar con la solución de ujo potencial. ParaRe = 2000 se tiene d = 0,0399.

7.5.2. Datos del problema

En cuanto a la resolución de la capa límite, una vez conocidas U1 y U2, se procede aconstruir el resto de condiciones de contorno. Teniendo en cuenta 7.1, se construye el primer

78 ETSII-UPM


perl de velocidades que utilizará el programa de resolución de las ecuaciones de la capalímite laminar como condición de contorno. Para ello, es necesario determinar primeramenteel valor máximo del vector y. Como se ha hecho con los anteriores programas de la placa planay del cilindro, se han ido probando distintos valores y observando los perles de velocidadesque se obtienen con cada uno de ellos. Finalmente se ha elegido el vector y de 0 a 10 con unpaso 4y = 0,01.

Se recuerda que y′′ =√ReR1y, por lo que si y′′max = 10, entonces se están resolviendo las

ecuaciones de la capa límite hasta una altura determinada y desde la pared inferior del vasoque varía en función del número de Reynolds. Por ello, el primer perl de velocidades es unaparábola según la ecuación 7.1, de tal manera que el máximo de la velocidad del primer perlserá

u′max =y′′max√Re

(2− y′′max√Re

)

Por otro lado, el vector x que recorre la longitud de la estenosis se ha elegido con un

paso 4x = 0,001 y longitud 2 long, donde long =√

(L2)2 + (R1 −R2)2, teniendo en cuenta

la geometría de la gura 7.2.

Dicho esto se puede construir el vector de velocidades Up a lo largo de la longitudde la estenosis a una altura ymax, valiendo el primer valor umax = u′maxU1 y el último U2.Entre ellos se hace una variación lineal como sugiere el artículo de Reese [10], todo elloadimensionalizando con R1y U1:

Up(x′) =

U2

U1− u′max2 longR1

x′ + u′max

Se ha contemplado la opción de construir el vector Up por conservación de caudales, esdecir,

Up(x) =U1R

21

2R(x)2

Esto no se puede realizar porque se recuerda que para cada caso se resuelve hasta unaaltura con dimensiones ymax distinta, y ésta se conserva a lo largo de la estenosis, es decir,en la garganta también se está resolviendo hasta la altura ymax, por lo que el caudal no seconserva según esta última ecuación. Consecuentemente, la variación lineal se ha dado porválida.

Además esta velocidad hay que proyectarla en dirección perpendicular a la pared delvaso para que sea correcta la resolución de las ecuaciones de la capa límite laminar, según semuestra en la gura 7.6:

Upcosα = UpL

2 long



α

α Up

Upcosα

Figura 7.6: Velocidad en los nodos superiores

El vector de presiones a lo largo de la estenosis se ha construido teniendo en cuenta laecuación de Bernouilli:

dp

dx= −Up

dUpdx

La validez de esta decisión radica en que los esfuerzos viscosos son importantes cercade las paredes, y fuera de la capa límite se pueden despreciar. Por ello la aproximación deBernouilli es prácticamente válida en la parte central de la geometría. Esta aproximacióntambién se realiza en el artículo de Reese.

Por lo tanto las condiciones de contorno para la velocidad horizontal U son:

u1(y) = u(x0, y) = u′p

UnN = Un

p

En cuanto a la velocidad vertical V , la condición de contorno es la misma que para laplaca plana y el cilindro:

vn0 = 0

7.5.3. Cálculo del esfuerzo cortante

Como ya se ha comentado, la comparación de los resultados se realizará a través del es-fuerzo cortante τ . Se recuerda que τ = µ∂u

∂y|y=0, por lo que con la variables adimensionalizadas

queda:

τ = µU1

√Re

R1

∂′u′

∂′y′′|y′′=0 ' µU1

√Re

R1

Un2 − Un

1

4y

Como se predice experimentalmente [10, 9], τmax se produce antes de la garganta, lo cualtambién se cumple con esta metodología, y se ve representando τ a lo largo de la estenosis.

80 ETSII-UPM


7.5.4. Resultados

Los resultados obtenidos se presentan a continuación grácamente en la gura 7.7,donde se pueden comparar con los experimentales. Además se presenta una comparación conlos obtenidos por Reese en su artículo en el cuadro 7.4.

102 103

Re

100

101

102

103

τ e

n di

nas/

cm2




Experimental

102 103

Re

100

101

102

103

τ e

n di

nas/

cm2





Experimental

Figura 7.7: Comparación de τmax calculado con los experimentales

Como puede apreciarse los resultados son muy similares a los obtenidos experimental-mente, lo cual da validez a este método. También se ha dado como resultado los puntos dedesprendimiento de la capa límite, que se presentarán más adelante.

Re M1 M2 M3 M4 M5 Fuente

1008,3 13,2 36,7 28,6 49,9 Experimental9,7 14,6 36,1 26,7 43,4 Calculado8,9 13,1 33,2 28,5 40,2 Reese



Cuadro 7.4: Comparación con los resultados experimentales y con los del artículo Reese(dina cm−2)

7.5.5. Comprobación de hipótesis

En esta sección se van a comprobar dos hipótesis tomadas y que están relacionadas.



En primer lugar, se procede a comprobar que el caudal de entrada al vasoQ1 = πR21U1

2=

π(R2 − dR2)2U2 es igual el caudal en la garganta, calculado como:

Q2 =

ˆ R2

0

u(r) 2πr dr

donde u(r) es la velocidad en la garganta. Teniendo en cuenta el cambio de variabley = R2 − r y poniendo las variables adimensionales queda:

Q2 =

ˆ 0

R′′2

2π(R2 − y′′R1√Re

)U1u′(y′′)

R1√Re

(−dy′′)

Retomando la notación sin comillas para las variables adimensionalizadas, se calculapor la regla de cuadratura del trapecio (5.2):

Q2 'R1√Re

2πU1

N∗−1∑i=1

4y2

(R2(u(i) + u(i+ 1))− R1√Re

(u(i) y(i) + u(i+ 1) y(i+ 1)))

donde N∗ es la posición del nodo que se encuentra en el eje del vaso sanguíneo, es decir,se calcula hasta R2 la integral. En cuanto a esto cabe comentar que hay casos en los queel vector de velocidades en la garganta está calculado hasta una altura menor que R2. Paraestos casos se ha creado un nuevo vector rellenando las componentes que faltan con velocidadU2 proyectada, pues se recuerda que el perl es uniforme en la garganta.

Por lo tanto, se ha comprobado que los caudales Q1 y Q2 son similares, como se muestraen la tabla 7.5. El error se ha calculado como |Q1−Q2|

Q1.

Como se puede ver los errores son aceptables teniendo en cuenta el error que conllevala aproximación de la integral por reglas de cuadratura y la resolución de las ecuaciones dela capa límite laminar por métodos numéricos, así como la suposición de la variación de lapresión según la ecuación de Bernouilli.

La segunda hipótesis a comprobar está relacionada con la anterior, y la vericación deque efectivamente el espesor de desplazamiento supuesto es igual al que se puede calcular enel perl, es decir

dR2 = δ∗ 'ˆ δ

0

(1− u

U) dy

donde δ se ha considerado conde uU2cosα

= 0,999.

De nuevo, teniendo en cuenta las adimensionalizaciones de variables, esta integral seha calculado por la regla de cuadratura del trapecio:

δ∗ 'ˆ δ

0

(1− u′U1

U2cosα)R1√Re

dy′′ ' R1√Re

N∗−1∑i=1

4y2

(1 +u′(i)U2

U1cosα

+ 1 +u′(i+ 1)U2

U1cosα

)

82 ETSII-UPM


Caso Re Q1(mLmin−1) Q2(mLmin

−1) Error (%)

M1100 115.45 105.72 8.4500 577.27 585.48 1.41000 1154.54 1127.89 2.3

M2100 115.45 116.96 1.3500 577.27 599.17 3.81000 1154.54 1135.58 1.6

M3100 115.45 126.01 9.1500 577.27 606.61 5.11000 1154.54 1134.11 1.8

M4100 115.45 116.47 0.9500 577.27 591.10 2.41000 1154.54 1116.38 3.3

M5100 115.45 123.91 7.3500 577.27 576.16 0.21000 1154.54 1080.58 6.4

Cuadro 7.5: Caudales

donde en este caso N∗ es la posición del nodo donde uU2cosα

= 0,999.

Los resultados obtenidos mostrados en la tabla 7.6 no son tan exactos como los de latabla 7.5. Sin embargo en todos lo casos la diferencia entre dR2 y δ∗es del mismo orden demagnitud y en ningún caso una es mayor que el doble que la otra. Como este resultado estárelacionado con el del caudal, pues ambos corroboran la validez del espesor de desplazamientosupuesto inicialmente, se consideran aceptables.

7.6. Extrapolación a mayores números de Reynolds: in-

troducción de la teoría del ujo potencial

La resolución de este problema se ha llevado a cabo de la misma manera que losproblemas de la placa plana y el cilindro. En primer lugar se ha resuelto el ujo potencial enFreeFem++ y posteriormente las ecuaciones de la capa límite laminar en Matlab.


El primer parámetro a determinar para este apartado es la velocidad de entrada alvaso sanguíneo del ujo. Para que luego sea congruente la comparación con el otro método,se ha elegido una velocidad tal que el perl de entrada es uniforme y de valor U1

2. Como se ha

explicado anteriormente, según la ecuación Re = ρU1R1

µ, esta velocidad depende del número de

Reynolds. Sin embargo, se recuerda que la resolución del problema en FreeFem++ requiere la



Caso Re δ∗(mm) dR2(mm)

M1100 0.631 0,471500 0.248 0,2621000 0.158 0,131

M2100 0.432 0.418500 0.187 0.2321000 0.127 0.116

M3100 0.243 0.315500 0.115 0.1751000 0.083 0.088

M4100 0.331 0.315500 0.155 0.1751000 0.112 0.088

M5100 0.186 0.315500 0.091 0.1751000 0.065 0.088

Cuadro 7.6: Espesor de desplazamiento

introducción de los datos adimensionalizados, por lo que la velocidad de entrada del ujo, quese adimensionaliza con ella misma por ser una de las variables elegidas para adimensionalizar,será U1 = 0,5.

Se ha resuelto la mitad del vaso sanguíneo, según se muestra en la gura 7.8. Conse-cuentemente, el eje del vaso sanguíneo se ha tomado como linea de corriente, de igual maneraque la pared del mismo. En esta gura aparece el valor de las funciones de corriente en cadafrontera así como la geometría considerada.

Figura 7.8: Geometría del problema

84 ETSII-UPM


El valor de las funciones de corriente es, por lo tanto, el siguiente:

ϕC1 = ϕC2 = ϕC3 = ϕC4 = 0

ϕC5 = ϕC7 = U1y

ϕC6 = U1R1

Cabe comentar que las presiones se han medido a una distancia h medida perpendicu-larmente a la pared del vaso, como ya se hizo en el programa del cilindro. Se tiene que cumplirque h sea mayor que el espesor de la capa límite δ. El problema es que no se puede conocerel espesor δ sin resolver las ecuaciones de la capa límite, para lo cual hace falta el vectorde presiones. Por lo tanto, se ha realizado un proceso iterativo de elección de h, resoluciónde ujo potencial y de las ecuaciones de la capa límite, y medida de δ en la garganta de laestenosis, en el perl de velocidad correspondiente a ese punto. Finalmente se ha llegado aque h = 0,12R2 y se ha comprobado para los cinco casos resueltos que h > δ. Como se verámás adelante en los perles de velocidades se puede medir el espesor de la capa límite adi-mensionalizado δ′′, que con dimensiones es δ. Se ha considerado δ′′ el nodo donde u

UN= 0,999.

En el cuadro 7.7 se muestran los valores de δ en la garganta para los distintos casos.

Caso h = 0,12R2 (mm) δthroat(mm)

M1 0,312 0,248M2 0,276 0,214M3 0,216 0,150M4 0,216 0,220M5 0,126 0,125

Cuadro 7.7: Espesor de la capa límite en la garganta

7.6.2. Capa límite

Para la resolución de las ecuaciones de la capa límite laminar por métodos numéricosse han introducido como datos el vector de presiones de ujo potencial.

En número de elementos del vector x es el mismo que el del vector de presiones queda FreeFem++, y a partir de eso se obtiene el valor de 4x. De esta manera, el número decomponentes ambos vectores es Nx. La longitud del vector de presiones es:

Nx =L

0,01

donde L es la longitud total de la estenosis, que en el programa se ha denominadolong, y 0,01 es el intervalo de medida de presión en el programa de FreeFem++. Si se



adimensionaliza con R1, la longitud adimensionalizada en el programa de Matlab será LR1,

por lo que igualando se tiene:

Nx =L

0,01=

L

R14x

Por lo tanto el paso del vector x es 4x = LR1Nx

, y variará en función del caso que seesté resolviendo.

Por otro lado, de igual manera que en caso de condiciones de contorno de Poiseuille, seha elegido el vector y de 0 a 10 con un paso de 4y = 0,01.

En cuanto a las condiciones de contorno, el primer perl de velocidades se ha tomadouniforme y de valor Um = U1

2, pues se está resolviendo con la condición de contorno de ujo

potencial y si se pusiera perl de velocidades de Poiseuille se estaría introducción el factorde la viscosidad en el problema. Además se ha puesto la condición de contorno 4.5.

La resolución de este problema ha dado como resultado los perles de velocidades, asícomo los puntos de desprendimiento de la capa límite.

7.7. Vericación del método

La validez de la aplicación de método con condiciones de contorno de ujo potencialradica en la comparación de los resultados con los obtenidos por el de condiciones de contornode ujo de Poiseuille, todo ello para números de Reynolds elevados. Como ya se ha explicado,el número de Reynolds al que se ha realizado es 2000 por ser el mayor posible que se incluyeen el régimen laminar. La comparación se lleva a cabo a través de la comparación de lospuntos de desprendimiento en los cinco casos estudiados. Para ello se presentan los puntosde desprendimiento obtenidos para Re = 1000 por el método de condiciones de contorno deujo de Poiseuille, que están contrastados con los datos experimentales, así como los puntosde desprendimiento por este mismo método para Re = 2000. Este último se ha resuelto deigual manera que los de Re = 100, 500, 1000 si bien en este caso d = 0,0399 como ya secomentó en la sección 7.5.1. Además se van a comparar con los puntos de desprendimientoobtenidos por el método de condiciones de contorno de ujo potencial.

Se espera que los puntos de desprendimiento del segundo método estén más adelantadosque los del primer método, pues el ujo potencial desprecia los efectos de la viscosidad, quees la responsable de la adherencia de la capa límite a la pared del vaso sanguíneo.

Por lo tanto, se comprueba que en todos los casos el punto de desprendimiento mástemprano es el de condiciones de contorno de la teoría potencial, seguido por el que ocurrecon Re = 2000 con condiciones de contorno de la teoría de Poiseuille. El más tardío es elque ocurre con Re = 1000 con la teoría de Poiseuille como condición de contorno. Además,teniendo en cuenta que las guras 7.9 y 7.10 tienen medidas reales en mm, se puede armarque para una misma longitud de estenosis (casos M1 y M2), el desprendimiento de la capa

86 ETSII-UPM


0 2 4 6 8 10 12 14L

0

0.5

1

1.5

2

2.5

3

3.5R

1

Puntos de desprendimiento para los casos M1 y M2

Pared del vaso M1

Pared del vaso M2

Flujo de Poiseuille

Flujo potencial

Flujo de Poiseuille con Re = 2000

0 5 10 15 20 25 30L

0

0.5

1

1.5

2

2.5

3

3.5

R1

Puntos de desprendimiento para los casos M3, M4 y M5

Pared del vaso M3

Pared del vaso M4

Pared del vaso M5

Flujo de Poiseuille

Flujo potencial


Figura 7.9: Puntos de desprendimiento de la capa límite

límite es más tardío para reducciones de la luz de la arteria menores, como es de esperar. Encuanto a igualdad de contracción del vaso (casos M3, M4 y M5), el desprendimiento es mástardío para el caso de menor longitud de la enfermedad (caso M5), siendo el más tempranode todos el de mayor longitud de estenosis (caso M4) evidentemente.

6.8 7 7.2 7.4 7.6 7.8 8 8.2 8.4 8.6L

0.7

0.8

0.9

1

1.1

1.2

1.3

R1

Puntos de desprendimiento para los casos M1 y M2

Pared del vaso M1

Pared del vaso M2

Flujo de Poiseuille

Flujo potencial


4 6 8 10 12 14L

1.45

1.5

1.55

1.6

1.65

1.7

1.75

R1

Puntos de desprendimiento para los casos M3, M4 y M5

Pared del vaso M3

Pared del vaso M4

Pared del vaso M5

Flujo de Poiseuille

Flujo potencial


Figura 7.10: Ampliaciones de los puntos de desprendimiento

7.8. Conclusiones

Se ha realizado el análisis de un tramo de una arteria de radio R1 = 3,5mm, en funcióndel grado de estenosis, analizando es valor de los esfuerzos cortantes. Un alto valor del mismopuede provocar el crecimiento de la placa de ateroma, pudiendo desprenderse y formar uncoágulo. La utilidad del análisis realizado radica en la posibilidad de predicción de situaciones



que no se puede observar experimentalmente, y como herramienta de ayuda a la comprensiónde esta patología que puede afectar al sistema cardiovascular.

Los dos métodos presentados se pueden emplear en el análisis del estado de la enfer-medad sanguínea, siendo más conservador el método que impone condiciones de contornode ujo potencial, si bien menos realista para esta aplicación. Se ha visto cómo afecta lageometría a los distintos puntos de desprendimiento de la capa límite, adelantándose conel crecimiento de la enfermedad. Se podría plantear un estudio paramétrico de la misma, sibien, observando las tendencias obtenidas, se llegaría a que no se produciría desprendimientopara L = 0 y H = 0, lo cual se reduce al absurdo de inexistencia de la enfermedad y vasoperfectamente cilíndrico.

Se recuerda de nuevo que tras el desprendimiento de la capa límite, ésta se reenganchaa la pared del vaso sanguíneo tras recorrer una cierta longitud, y que durante ésta, se puedenoriginar recirculaciones de provoquen la deposición de plaquetas en la pared del vaso, crecien-do la placa de ateroma. En cuanto a esto, es importante tener en cuenta que los dos métodospresentados nos son válidos para el estudio de ujo sanguíneo posterior al desprendimientode la capa límite.

También se ha comentado la importancia de la aparición de la disfunción endotelial siel esfuerzo cortante máximo supera el valor de 380 dina cm−2, lo cual ocurre para númerosde Reynolds altos (en torno a Re = 1000) y reducciones de la luz de la arteria de aproxi-madamente el 75 %, todo ello según el cuadro 7.2. Todos estos casos se incluyen en estenosissevera-grave según la OMS.

7.9. Mejoras en la geometría

La geometría de la estenosis real evidentemente no es la que se ha supuesto en laresolución del problema. Esta patología se diagnostica a partir de angiografías, y en la gura7.1 se muestra una imagen de una angiografía de una arteria coronaria donde se apreciaclaramente el estrechamiento del vaso. Como puede verse la pared de la arteria es muyirregular en la zona de la estenosis.

La geometría utilizada es la de la gura 7.2 si bien no es realista. Existen otras funcionesmatemáticas que aproximan mejor la forma de la estenosis como pueda ser [7]:

r(x) = R1 −H

2[1 + cos(π(1− 2x

L))]

Se ha representado esta función en la gura 7.11 para el caso M4, y se puede ver quees mucho más cercana a la geometría que se presenta realmente en los vasos afectados. Estapodría ser una de las lineas de trabajo futuras de este proyecto para obtener resultados másrealistas sobre esta enfermedad.

88 ETSII-UPM


0 5 10 15 20 25

x

-3

-2

-1

0

1

2

3

r(x)

Figura 7.11: Geometría de la estenosis más realista



90 ETSII-UPM

Conclusiones generales

Como se ha visto a lo largo de todo este proyecto, el tratamiento de las ecuaciones querigen la Mecánica de Fluidos no es sencillo. Su resolución es una herramienta muy útil paraconocer el comportamiento de los uidos. Para ello ha sido necesaria una previa simplicaciónde las ecuaciones dividiendo el ujo en dos regiones bien diferenciadas en función de a laimportancia de los efectos viscosos en las mismas.

Las teorías de ujo potencial y de capa límite laminar son la base de este proyecto.Ambas teorías son muy poderosas ya que simplican las ecuaciones de Navier Stokes parapoder resolver el campo uido.

Sin embargo, y como es evidente, tienen limitaciones, especialmente en el campo deaplicación, pues la teoría de ujo potencial, al despreciar los efectos viscosos requiere númerosde Reynolds del ujo muy altos. Por otro lado, la teoría de la capa límite laminar, como supropio nombre indica, se aplica a régimen laminar, lo que supone una limitación a la resolucióndel régimen de ujo turbulento. Otra de las limitaciones de la teoría de la capa límite laminares que no permite el estudio del ujo a partir del punto de desprendimiento de la misma.

El tratamiento de ambas teorías se ha realizado gracias a los métodos numéricos. Comoes sabido, su aplicación lleva asociada un error de cálculo determinado. Éste depende en granparte del modelo de aproximación de las derivadas parciales que aparecen en las ecuaciones deNavier Stokes, así como del tamaño de malla. Para la resolución de las ecuaciones de la capalímite laminar el modelo en derivadas parciales y el mallado se han elegido personalmenteteniendo en cuenta modelos propuestos para la resolución de estas ecuaciones [12].

Otro de los aspectos más importantes de este proyecto es que en todo momento sehan ido comparando los resultados obtenidos con resultados bien analíticos o bien experi-mentales consultados en la bibliografía. Esto ha permitido el control del error cometido entodo momento y consecuentemente la mejora del modelo matemático de manera constante yprogresiva a lo largo del proyecto

Cabe resaltar que las ecuaciones de Navier Stokes se pueden resolver directamentea través de códigos Dinámica de Fluidos Computacional (CFD). Sin embargo, una de lasventajas del modelo matemático que concilia las teorías de ujo potencial y capa límitelaminar es su simplicidad tanto en planteamiento como en tratamiento. A lo largo de todoel proyecto ha sido posible relacionar la fenomenología física asociada a cada situación con

91


el problema matemático que lo respalda. Un claro ejemplo de esto es la justicación de laelección de unas condiciones de contorno para cada caso y geometría o la elección de valorde las funciones de corriente.

En cuanto al desarrollo de la aplicación, se concluye y se demuestra la validez del modelomatemático desarrollado para el estudio del punto de desprendimiento de la capa límite, elcual genera disipación de energía y no es deseable. Como se ha visto, unas ecuaciones tansimples como aquellas que se han manejado permiten un amplio estudio de posibilidades(geometría, tipos de ujo,...) todo ello a pesar de las limitaciones del método ya comentadas.Además el tiempo de cálculo de los programas es muy inferior al que requiere la resolucióncompleta de las ecuaciones de Navier Stokes con CFD. De hecho el tiempo de cálculo es tresminutos para cada programa realizado frente a 12 horas que tardarían los códigos CFD.

En cuanto a posibles lineas futuras del proyecto, se contemplan varias opciones:

Mejora del modelo matemático por renamiento de la malla así como del sistema deecuaciones de diferencias nitas a resolver en la teoría de capa límite laminar. De estamanera se disminuiría el error cometido, aumentando la precisión en el cálculo.

Comparación de los resultados obtenidos con los que resultan de la resolución de lasecuaciones de Navier Stokes completas, para lo cual sería necesario el empleo de pro-gramas como FreeFem++ o OpenFOAM.

Desarrollo en profundidad de la succión y de otros métodos de control de la capa límitecomo el soplado, inyección de un uido o la alteración de la supercie del cuerpo.

Aplicaciones variadas que se enmarquen en el rango de números de Reynolds del régimenlaminar. En relación a la aplicación en medicina de la sangre, se contemplaría la opciónde resolver el caso de ujo alrededor de la aneurisma arterial, que consiste en unadilatación del vaso sanguíneo en una determinada zona.

Ampliación a capa límite turbulenta.

Como es lógico, existen más modelos que permiten resolver el ujo alrededor de un cuerpoademás del que ha desarrollado en este proyecto. Sin embargo, la sencillez de su planteamientoy el reducido tiempo de cálculo lo hacen idóneo para la compresión de muchos fenómenosasociados a los uidos. A lo largo de todo el proyecto se ha podido entender la relación entrela física del problema y la matemática utilizada para su resolución.

92 ETSII-UPM

Apéndice A

Planicación temporal y presupuesto

A.1. Planicación temporal

Este apartado tiene como objetivo mostrar el conjunto de actividades que se han llevadoa cabo durante la realización de este Trabajo de Fin de Grado, así como su duración y relación.Para ello se presentan a continuación el diagrama de Gantt y la Estructura de Descomposicióndel Proyecto (EDP).

A.1.1. Estructura de Descomposición del Proyecto

La EDP recoge las tareas a realizar de forma estructurada y disciplinada. Para ellose tienen en cuenta todas las etapas del proyecto, desde los niveles superiores a los másinferiores, llamados paquetes de trabajo. Para el caso de este proyecto, la EDP se muestraen la gura A.1.

Como puede observarse, se ha dividido en cuatro ramas diferenciadas que engloban lossiguientes campos:

Formación previa: como es evidente, incluye la familiarización con los conceptos básicosdel proyecto como puedan ser las teorías empleadas en el mismo o las ecuaciones deNavier Stokes. Además se incluye el curso realizado para aprender a usa el softwareFreeFem++.

Desarrollo del modelo matemático: incluye campos de trabajo como la manipulaciónde las ecuaciones de acuerdo con los modelos a realizar, así como el desarrollo e imple-mentación del código fuente en ambos programas. Además se realizaron una serie desimulaciones de prueba que dan pie la siguiente rama.

93

APÉNDICE A. PLANIFICACIÓN TEMPORAL Y PRESUPUESTO

Vericación del modelo matemático: para ello se han probado dos geometrías y se hancontrastado con los resultados que aparecen en la literatura.

Aplicación: incluye fases como la documentación y el desarrollo de los programas per-tinentes para llevar a cabo la simulación, así como la obtención de resultados.

A.1.2. Diagrama de Gantt

Se trata de una herramienta que muestra la duración en el tiempo dedicado a cadaactividad así como la relación existente entre ellas. Representa de forma gráca el progresodel proyecto, e incluye relaciones del tipo:

Fin a comienzo: aquellas actividades que para comenzar requieren de la nalizaciónde otras actividades previas. Un claro ejemplo es la programación en Matlab y Free-Fem++, pues ha de realizarse cuando ya se haya recopilado información suciente sobreel software y las ecuaciones uidodinámicas.

Comienzo a comienzo: el comienzo de una actividad está marcado por el inicio de otra.Por ejemplo la vericación de resultados comienza una vez que se han obtenido losmismos.

Se muestra en la gura A.2.

A.2. Presupuesto del proyecto

Como es previsible, la ejecución de este proyecto no requiere de inversión alguna, puesesencialmente se han llevado a cabo simulaciones de tipo teórico y no hay coste material realen sí mismo. Sin embargo, sí es posible estimar el coste total del proyecto teniendo en cuentadistintos puntos:

Coste de personal

Incluye las horas personales dedicadas al proyecto. Además se podría incluir el coste delasesoramiento por parte de los tutores a través del precio por hora de los mismos, que seestimaría en 80 ¿/h.

Horas trabajadas Precio por hora (¿/h) Precio total (¿)

380 40 13600Total 13600 ¿

94 ETSII-UPM


Coste de formación teórica

En este apartado se incluye el coste de la asistencia al curso de introducción a FreeFem++impartido en la asignatura Ingeniería de Fluidos del Máster en Ingeniería Industrial. Se tratade una asignatura de 4.5 créditos, y el precio por crédito es 62 ¿. Se ha asistido a 18 horas delas 42 que incluye la asignatura, por lo que el precio es el que aparece en la siguiente tabla:

Horas de asistencia al curso Precio por hora de curso Precio (¿)

18 6.64 119.6Total 119.6 ¿

Coste de herramientas informáticas

Este punto incluye el coste de los equipos informáticos y programas, teniendo en cuenta suamortización y el tiempo de uso.

Producto Precio (¿) Amortización Tiempo de uso Precio total (¿)

Ordenador ASUS 319 4 años 9 meses 59.81Licencia Matlab 69 1 año 9 meses 12.94

Licencia FreeFem++ 0 - 9 meses 0Total 72.75 ¿

Coste de material

Se incluye el coste de los principales artículos utilizados para realizar la aplicación teórica delproyecto.

Producto Número Precio (¿)

Artículos 4 30Total 30 ¿

Por lo tanto, se tiene nalmente el presupuesto del proyecto:

Concepto Precio (¿)

Coste de personal 13600Coste de formación teórica 119.6

Coste de herramientas informáticas 72.75Coste de material 30

TOTAL 13822.32 ¿

Cuadro A.1: Presupuesto del proyecto



Figura A.1: EDP

96 ETSII-UPM


Figura A.2: Diagrama de Gantt



98 ETSII-UPM

Apéndice B

Reexión general sobre aspectos deresponsabilidad social del trabajo

A lo largo de este trabajo se han utilizado diversas herramientas de ingeniería con elobjetivo de simular una situación que se atribuye a la medicina de la sangre. Hoy en día lassimulaciones tienen una gran impacto e importancia en el mundo de la ciencia. La ingenieríaque se está desarrollando se basa en simulaciones que abarcan desde la operación de unreactor nuclear o el vuelo de un avión hasta el funcionamiento del corazón. Concretamente,la Mecánica de Fluidos ha ganado terreno en el mundo de las simulaciones gracias, en parte,al reciente desarrollo de programas de Dinámica de Fluidos Computacional (CFD). Por otrolado, los recientes avances médicos son fruto de muchos años de experimentación y simulación,por lo que se puede decir que el avance de la tecnología permite la mejora de las condicionesde vida.

Cabe comentar que según la Sociedad Española de Cardiología (SEC), las enfermedadescardiovasculares provocan 65 veces más fallecimientos que los accidentes de tráco. Ademásla SEC arma que el 80% de estas enfermedades se podrían prevenir con un estilo de vidasaludable. Es por ello que la dedicación de recursos al estudio de este tipo de patologías estájusticada, y la ingeniería de hoy en día está orientada hacia una investigación de caráctermás ético que las que se realizaban hace un siglo.

La aplicación de las ecuaciones de la Mecánica de Fluidos a la medicina de la sangrepermite el análisis de muy variadas situaciones como pueda ser los aspectos hidrodinámicosde las válvulas utilizadas en cirugía cardíaca, los fenómenos de ondas de presión, las recircu-laciones, el ujo turbulento o la interacción entre la sangre y dispositivos que se introducenen el sistema circulatorio.

Una de las consecuencias más inmediatas del amplio uso de las simulaciones es el ahorroen material que supone que el hecho de poder analizar una determinada situación. En esteproyecto la simulación de la geometría de la estenosis ha permitido analizar la gravedad de laenfermedad en base a la contracción del vaso sanguíneo y a la magnitud del esfuerzo cortante.

99

APÉNDICE B. REFLEXIÓN GENERAL SOBRE ASPECTOS DE RESPONSABILIDADSOCIAL DEL TRABAJO

Este parámetro uidodinámico es difícil de medir en un paciente real con estenosis arterial,y sin embargo, los expertos e investigadores han conado en la mecánica computacional y losmodelos analíticos para estimar su valor.

Además las simulaciones permiten variar las condiciones en que se produce un fenó-meno, abriéndose un amplio camino hacia la experimentación. Para el caso de la estenosis,basándose en datos experimentales y casos conocidos, ha sido posible variar la geometría dela arteria o incluso el carácter del ujo. Como se comentó en el Capítulo VII, el ujo pulsátilsería más real que el estacionario, y se podría aproximar la geometría de la estenosis a unamás realista a través de modelos matemáticos.

En cuanto a los aspectos sociales y éticos, está claro que el impacto de las simulacionesen el campo de la medicina es muy grande, si bien hay que tener en cuenta, como es evidente,que la incertidumbre existente en lo que concierne al cuerpo humano es alta. En este sentido,la responsabilidad social y ética es importante ya que se ha de ser consciente de que cualquierdesviación del modelo teórico o de lo previsto puede ser fruto de una desviación en el cursode la enfermedad motivada por factores externos a la misma, y no de un error del modelomatemático.

Es por todo ello que los estudios de bioingeniería son un apoyo a la medicina, unaherramienta para su avance y mejora, que no deja de crecer.

100 ETSII-UPM

Apéndice C

Código fuente

C.1. Placa plana

Se presentan los códigos de Matlab, ya que los de FreeFem++ se presentarán másadelante para otras geometrías.

C.1.1. Solución pseudoanalítica

1 %−−−−−−−Función para poner la ecuación diferencial en forma explícita según2 %criterio de Runge−Kutta−−−−−−−%3 function df=placa(~,f)4 df=zeros(3,1);5 df(1)=f(2);6 df(2)=f(3);7 df(3)=−f(1)*f(3)/2;8 end

1 %−−−−−−−−−−−−−−−−−−−Programa método del disparo−−−−−−−−−−−−−−−−−−−%2 %Intervalo de eta3 interveta=0:0.01:10;4 %Vectores para almacenar datos5 q=zeros(100);6 r=zeros(100);7 %Bucle para hallar el valor de la condicion inicial de f'' a partir de la8 %condición conocida de que f'(infinito)=19 for i=1:100

10 q(i)=i/100;11 %La condición de contorno sobre f'' irá de 0.01 a 112 [~,f]=ode45(@placa,interveta,[0 0 q(i)]);13 r(i)=f(length(f),2);

101

APÉNDICE C. CÓDIGO FUENTE

14 %En el vector r se guardan lo valores de f' a lo largo el bucle.15 end16 plot(q,r)17 xlabel('Condición de contorno')18 ylabel('$ \fracdfd\eta$','Interpreter','latex');19 grid on20 %En la gráfica se ven los valores de la condición f' en función de los de21 %f''. Se está buscando una condición sobre f'' de tal manera que se cumpla22 %f'=1. En la gráfica se ve que esto ocurre para q=0.332, r=1, es decir,23 %para f'=1 f''=0.332

1 %−−−−−Programa resolución ecuación pseudoanalítica. Perfil de Blasius−−−−−%2 function[u]=blasius(Up)3 etamax=7;4 interveta=0:0.01:etamax;5 %Solver ode456 [eta,f]=ode45(@placa,interveta,[0 0 0.332]);7 %Representación de la solución8 figure9 plot(eta,f(:,2));

10 xlabel('$\eta$','Interpreter','latex');11 ylabel('$ \fracdfd\eta$','Interpreter','latex');12 title('Curva teórica de la distribución de velocidades en la capa límite

laminar');13 grid on14 axis([0,5,0,1]);15

16 u=Up*f(:,2);17 y=eta;18 figure19 plot(u,y)20 xlabel('U')21 ylabel('Altura')22 title('Perfil de velocidades en la capa límite para una placa plana')23 grid on24

25 %Guardo para comparar porteriormente con lo obtenido por métodos numéricos26 fprima=f(:,2);27 save('f_eta','fprima','eta','y','u')28 end

C.1.2. Solución por métodos numéricos

1 %−−−−−−−Programa resolución capa límite en placa plana sin succión−−−−−−−%2 function[U,V]=mn_placa(Re,Up)3 Dx=0.0001;4 Dy=0.01;5 x=0:Dx:1;

102 ETSII-UPM


6 y=0:Dy:7;7 N=length(y);8 Nx=length(x);9 U=zeros(N,Nx);

10 V=zeros(N,Nx);11 uold=ones(N,1); %Primer perfil de velocidades recto y de valor 112 uold=Up*uold;13 vold=zeros(N,1);14 U(:,1)=uold;15 V(:,1)=vold;16 Cf=zeros(Nx,1);17 Cfteor=zeros(Nx,1);18 A=zeros(N,N);19 B=zeros(N,1);20 K=zeros(N,N);21 L=zeros(N,1);22

23 %Bucle que recorre la placa24 for n=2:Nx25 %Construcción de las matrices A y B26 for i=2:N−127 A(i,i−1)=−vold(i)/Dy−1/Dy^2;28 A(i,i)=vold(i)/Dy+uold(i)/Dx+2/Dy^2;29 A(i,i+1)=−1/Dy^2;30 B(i,1)=(uold(i))^2/Dx;31 end32 %Condiciones de contorno33 A(1,:)=0; A(1,1)=1;34 A(N,:)=0; A(N,N)=1;35 B(1,1)=0; B(N,1)=1;36 %Resolución del primer sistema37 Usol=A\B;38 Usol(1)=0;39 U(:,n)=Usol;40 %Construcción de las matrices K y L41 for i=2:N42 K(i,i−1)=−1/Dy;43 K(i,i)=1/Dy;44 L(i,1)=−(U(i,n)−uold(i,1))/Dx;45 end46 %Condiciones de contorno47 K(1,:)=0; K(1,1)=1;48 L(1)=0;49 %Resolución del segundo sistema50 Vsol=K\L;51 Vsol(1)=0;52 V(:,n)=Vsol;53 %Se guardan las variables del nodo para utilizarlas en el siquiente54 %bucle55 uold=Usol;56 vold=Vsol;57 end58

59 %Representacion de U



60 figure61 plot(x(1)+U(:,1),y);62 title('Perfiles de velocidad horizontal en diferentes puntos de la placa');63 xlabel('U');64 ylabel('Altura');65 hold on66 grid on67 for i=1:n68 if mod(i,200)==069 plot(x(i)+U(:,i),y)70 end71 end72 plot(x(Nx)+U(:,Nx),y);73

74 %Representacion de V75 figure76 plot(x(1)+V(:,1),y);77 title('Perfiles de velocidad vertical en diferentes puntos de la placa');78 xlabel('V');79 ylabel('Altura');80 hold on81 grid on82 for i=1:n83 if mod(i,200)==084 plot(x(i)+V(:,i),y)85 end86 end87 plot(x(n)+V(:,n),y)88

89 %Superposicion con la solución con la pseudoanalítica. Gráfica 190 upseudo=load('f_eta.mat');91 figure92 plot(upseudo.eta,upseudo.fprima,'b');93 xlabel('$\eta$','Interpreter','Latex')94 ylabel('$ \fracdfd\eta$','Interpreter','latex');95 grid on96 hold on97 eta=y;98 plot(eta,U(:,Nx)/Up,'r');99 axis([0,5,0,1]);

100 legend('Solución pseudoanalítica','Solución numérica','Location','southeast')

101

102 %Superposición con la solución con la pseudoanalítica. Gráfica 2103 upseudo=load('f_eta.mat');104 figure105 plot(upseudo.u,upseudo.y,'b');106 xlabel('U'); ylabel('Altura');107 grid on108 hold on109 plot(U(:,Nx),y,'r');110 legend('Solución pseudoanalítica','Solución numérica','Location','southeast

')111

104 ETSII-UPM


112 %Coeficiente de fricción113 for i=2:Nx114 Cf(i,1)=2/sqrt(Re)*(U(2,i)−U(1,i))/Dy;115 Cfteor(i,1)=0.664/sqrt(Re*x(i));116 end117

118 figure119 loglog(x(2:Nx),Cf(2:Nx),'r');120 hold on121 loglog(x(2:Nx),Cfteor(2:Nx),'b');122 title('Coeficiente de fricción en ejes logarítmicos');123 xlabel('Longitud adimensionalizada de la placa');124 ylabel('$logC_f(x)$','interpreter','latex');125 %La primera solucion no tiene sentido pintarla porque viene del infinito y126 %no hay fricción127 grid on128 legend('Solución pseudoanalítica','Solución numérica')129

130 figure131 plot(x(2:Nx),Cf(2:Nx),'r');132 hold on133 plot(x(2:Nx),Cfteor(2:Nx),'b');134 title('Coeficientes de fricción teórico y numérico');135 xlabel('Longitud adimensionalizada de la placa');136 ylabel('$C_f(x)$','interpreter','latex');137 legend('Solución pseudoanalítica','Solución numérica')138 grid on139

140 %Valor medio del coeficiente de fricción141 sum=0;142 for i=2:Nx143 sum=sum+Dx/2*(U(2,i)−U(1,i)+U(2,i−1)−U(1,i−1))/Dy;144 end145 Cfmed=2/sqrt(Re)*sum;146 fprintf('Coeficiente de fricción medio Cf*sqrt(Re)=%d\n',Cfmed*sqrt(Re))147

148 %Diferencia entre Cf teótico y Cf numético. Error149 error=zeros(Nx);150 for i=2:Nx151 error(i)=(Cfteor(i,1)−Cf(i,1))/Cfteor(i,1)*100;152 end153 figure154 plot(x(2:Nx),error(2:Nx));155 title('Error entre coeficientes de fricción teórico y numérico en %');156 grid on157

158 %Se guardan los datos para comparar con el caso con succión159 Usinsucc=U(:,Nx);160 save('usinsucc','Usinsucc')161 end

1 %−−−−−−−Programa resolución capa límite en placa plana con succión−−−−−−−%



2 function[U,V]=mn_placa_succion(Re,Up,t)3 p=load('presiones_succion.txt');4 %Se carga el vector de presiones de Free Fem5 Nx=length(p)−1;6 Dy=0.01;7 y=0:Dy:7;8 Dx=1/Nx;9 x=0:Dx:(1−Dx);

10 dp=diff(p)/Dx;11 vs=−t*Up; %Velocidad de succión real, sin adimensionalizar12 N=length(y);13 Nx=length(x);14 U=zeros(N,Nx);15 V=zeros(N,Nx);16 uold=ones(N,1);17 uold=Up*uold;18 vold=zeros(N,1);19 U(:,1)=uold;20 V(:,1)=vold;21 A=zeros(N,N);22 B=zeros(N,1);23 K=zeros(N,N);24 L=zeros(N,1);25

26 %Bucle que recorre la placa;27 for n=2:Nx28 %Construcción de las matrices A y B29 for i=2:N−130 A(i,i−1)=−vold(i)/Dy−1/Dy^2;31 A(i,i)=vold(i)/Dy+uold(i)/Dx+2/Dy^2;32 A(i,i+1)=−1/Dy^2;33 B(i,1)=−dp(n)+(uold(i))^2/Dx;34 end35 %Condiciones de contorno36 A(1,:)=0; A(1,1)=1;37 A(N,:)=0; A(N,N)=1;38 B(1,1)=0; B(N,1)=1;39 %Resolución del primer sistema40 Usol=A\B;41 Usol(1)=0;42 U(:,n)=Usol;43 %Construcción de las matrices K y L44 for i=2:N45 K(i,i−1)=−1/Dy;46 K(i,i)=1/Dy;47 L(i,1)=−(U(i,n)−uold(i,1))/Dx;48 end49 %Condiciones de contorno50 K(1,:)=0; K(1,1)=1;51 L(1)=vs*sqrt(Re)/Up; %Velocidad de succión adimensionalizada52 %Resolución del segundo sistema53 Vsol=K\L;54 Vsol(1)=0;55 V(:,n)=Vsol;

106 ETSII-UPM


56 %Se guardan las variables del nodo para utilizarlas en el siquiente57 %bucle58 uold=Usol;59 vold=Vsol;60 end61

62 %Representacion de U−−igual que en el programa sin succión63

64

65 %Evolución de los perfiles en distintos puntos de la placa66 figure67 plot(−vs*sqrt(Re)/Up*y,U(:,4)/Up)68 hold on69 plot(−vs*sqrt(Re)/Up*y,U(:,11)/Up)70 plot(−vs*sqrt(Re)/Up*y,U(:,42)/Up)71 plot(−vs*sqrt(Re)/Up*y,U(:,92)/Up)72 plot(−vs*sqrt(Re)/Up*y,U(:,163)/Up)73 plot(−vs*sqrt(Re)/Up*y,U(:,253)/Up)74 plot(−vs*sqrt(Re)/Up*y,U(:,486)/Up)75 grid on76 title('Perfil de velocidades en distintos puntos de la placa')77 xlabel('$\frac−y v_s\nu$','Interpreter','latex')78 ylabel('$\fracuUp$','Interpreter','latex')79 t=legend('$\sqrt\xi=0.1$','$\sqrt\xi=0.2$','$\sqrt\xi=0.4$','$\sqrt\xi=0.6$

','$\sqrt\xi=0.8$','$\sqrt\xi=1$','$\sqrt\xi=1.4$','Location','southeast');

80 set(t,'Interpreter','Latex')81 axis([0,3.6,0,1])82

83 %Comparación con los perfiles sin succión84 Usin=load('usinsucc.mat');85 figure86 plot(y/1.7208,Usin.Usinsucc/Up,'r')87 hold on88 plot(y,(1−exp(−y)),'g')89 plot(−vs*sqrt(Re)/Up*y,U(:,Nx)/Up,'b')90 grid on91 legend('Perfil sin succión','Perfil teórico con succión','Perfil obtenido

con succión en el final de la placa','Location','southeast')92 title('Comparación de los perfiles con y sin succión')93 xlabel('$\frac−y v_s\nu$','Interpreter','latex')94 ylabel('$\fracuUp$','Interpreter','latex')95 axis([0,4,0,1])96 end



C.2. Cilindro

C.2.1. Sin succión

1 //−−−−−−Programa resolución flujo potencial en cilindro sin succión−−−−−−//2 //Datos3 real r=3; //radio del cilindro4 real w=4*r; //altura del mallado5 real U=1;6

7 //Construcción del mallado8 mesh Th;9 border C1(t=0,1)x=w;y=3*w*t;;

10 border C2(t=1,−1)x=w*t;y=3*w;;11 border C3(t=1,0)x=−w;y=3*w*t;;12 border C4(t=−w,−r)x=t;y=0;;13 border C5(t=pi,0)x=r*cos(t);y=r*sin(t);;14 border C6(t=r,w)x=t;y=0;;15 Th=buildmesh(C1(100)+C2(50)+C3(100)+C4(50)+C5(pi*100*r)+C6(50));16 plot(Th,wait=0,cmm=" Mallado");17

18 //Resolución del problema19 //Definición del espacio de elementos finitos20 fespace Vh2(Th,P2);21 Vh2 phi,g;22 fespace Vh1(Th,P1);23 Vh2 vx,vy,p;24 //Resolución de la ecuación lapl(phi)=025 problem laplace(phi,g)=int2d(Th)(dx(phi)*dx(g))+int2d(Th)(dy(phi)*dy(g))26 +on(C4,C5,C6,phi=0)27 +on(C1,C3,phi=U*y)28 +on(C2,phi=U*3*w);29 laplace;30 //Representación de la solución31 plot(cmm=" Funcines de corriente ",phi,wait=0,nbiso=30);32 vx=dy(phi);33 vy=−dx(phi);34 plot(cmm=" Velocidades en x y en y",[vx,vy],wait=0);35 p =0.5*(1−vx^2 − vy^2);36 plot(cmm=" Campo de presiones ",p,ColorScheme=1,fill=1,nbiso=20,wait=0);37 //Lectura de la presión a una distancia h de la semicircunferencia38 real h=0.025*r;39 ofstream gp("presiones.txt");40 for (real i=0;i<(pi*r);i=i+0.01)41 42 gp<<p(−(r+h)*cos(i/r),(r+h)*sin(i/r))<<endl;43

108 ETSII-UPM


C.2.2. Con succión

1 //−−−−−−Programa resolución flujo potencial en cilindro sin succión−−−−−−//2 //Datos3 real r=3; //radio del cilindro4 real w=4*r;5 real U=1;6 real vs=0.02*U;7 real Us=U−vs*pi*r/(3*w);8

9 //Construcción del mallado10 mesh Th;11 border C1(t=0,1)x=w;y=3*w*t;;12 border C2(t=1,−1)x=w*t;y=3*w;;13 border C3(t=1,0)x=−w;y=3*w*t;;14 border C4(t=−w,−r)x=t;y=0;;15 border C5(t=pi,0)x=r*cos(t);y=r*sin(t);;16 border C6(t=r,w)x=t;y=0;;17 Th=buildmesh(C1(50)+C2(50)+C3(50)+C4(50)+C5(pi*50*r)+C6(50));18 plot(Th,wait=0,cmm=" Mallado");19

20 //Resolución del problema21 fespace Vh2(Th,P2);22 Vh2 phi,g;23 fespace Vh1(Th,P1);24 Vh2 vx,vy,p;25

26 problem laplace(phi,g)=int2d(Th)(dx(phi)*dx(g))+int2d(Th)(dy(phi)*dy(g))27 +on(C1,phi=Us*(y−3*w))28 +on(C2,phi=0)29 +on(C3,phi=U*(y−3*w))30 +on(C4,phi=−3*w*U)31 +on(C5,phi=−3*w*U+vs*r*atan(y,−x))32 +on(C6,phi=−3*Us*w);33 laplace;34 plot(cmm=" Funcines de corriente ",phi,wait=0,nbiso=100);35 vx=dy(phi);36 vy=−dx(phi);37 plot(cmm=" Velocidades en x y en y",[vx,vy],wait=0);38 p =0.5*(1−vx^2 − vy^2);39 plot(cmm=" Campo de presiones ",p,ColorScheme=1,fill=1,nbiso=20,wait=0);40

41 real h=0.025*r;42 ofstream gp("presiones_succion.txt");43 for (real i=0;i<(pi*r);i=i+0.01)44 45 gp<<p(−(r+h)*cos(i/r),(r+h)*sin(i/r))<<endl;46

El programa de resolución de las ecuaciones de la capa límite es el mismo, y se muestra acontinuación.



1 %−−−−−−−Programa resolución capa límite en cilindro con/sin succión−−−−−−−%2 function[U,V]=mn_cilindro(Re,Up,t)3 %Si t=1 con succión, si t=0 sin succión4 if t==05 vs=0;6 p=load('presiones.txt');7 else8 vs=−t*Up;9 p=load('presiones_succion.txt');

10 end11 Nx=length(p)−1;12 Dy=0.01;13 y=0:Dy:3.5;14 Dx=pi/Nx;15 x=0:Dx:(pi−Dx);16 dp=diff(p)/Dx;17

18 figure19 plot(x,dp)20 grid on21 title('Gradiente de presiones')22 xlabel('Longitud de la circunferencia adimensionalizada')23

24 N=length(y);25 U=zeros(N,Nx);26 V=zeros(N,Nx);27 uold=zeros(N,1);28 vold=zeros(N,1);29 U(:,1)=uold;30 V(:,1)=vold;31 A=zeros(N,N);32 B=zeros(N,1);33 K=zeros(N,N);34 L=zeros(N,1);35

36 %Bucle que recorre la semicircunferencia37 for n=2:Nx38 for i=2:N−139 A(i,i−1)=−vold(i)/Dy−1/Dy^2;40 A(i,i)=vold(i)/Dy+uold(i)/Dx+2/Dy^2;41 A(i,i+1)=−1/Dy^2;42 B(i,1)=−dp(n)+(uold(i))^2/Dx;43 end44 A(1,:)=0; A(1,1)=1;45 A(N,:)=0; A(N,N)=1;46 B(1,1)=0; B(N,1)=2*Up*sin(x(n));47 %Resolución del primer sistema48 Usol=A\B;49 Usol(1)=0;50 U(:,n)=Usol;51

52 for i=2:N53 K(i,i−1)=−1/Dy;

110 ETSII-UPM


54 K(i,i)=1/Dy;55 L(i,1)=−(U(i,n)−uold(i,1))/Dx;56 end57 K(1,:)=0; K(1,1)=1;58 L(1)=vs*sqrt(Re)/Up;59 %Resolución del segundo sistema60 Vsol=K\L;61 Vsol(1)=0;62 V(:,n)=Vsol;63 %Se guardan las variables del nodo para utilizarlas en el siquiente64 %bucle65 uold=Usol;66 vold=Vsol;67 %Se deja de resolver cuando se desprende la capa límite68 if U(10,n)<069 break70 end71 end72 if U(10,n)<073 fprintf('Para Re=%d \n',Re)74 fprintf('Se separa en phi=%d grados \n',x(n)*180/pi)75 else76 n=Nx;77 end78

79 %Representacion de U80 figure81 plot(x(1)+U(:,1),y);82 xlabel('U');83 title('Perfiles de velocidad horizontal en diferentes puntos de la

semcircunferencia');84 ylabel('Altura');85 hold on86 grid on87 for i=1:n88 if mod(i,20)==089 plot(x(i)+U(:,i),y)90 end91 end92 plot(x(n)+U(:,n),y)93

94 %Representacion de V95 figure96 plot(x(1)+V(:,1),y);97 xlabel('V');98 title('Perfiles de velocidad vertical en diferentes puntos de la

semicircunferencia');99 ylabel('Altura');

100 hold on101 grid on102 for i=1:n103 if mod(i,20)==0104 plot(x(i)+V(:,i),y)105 end



106 end107 plot(x(n)+V(:,n),y)108

109 %Comparación con las curvas Schlihting (solución analítica)110 figure111 xlabel('$\fracuU$','Interpreter','latex')112 ylabel('$\fracyR\sqrt\fracU R\nu$','Interpreter','latex');113 hold on114 grid on115 plot(U(:,floor(1+pi*20/(Dx*180)))/(2*Up*sin(20*pi/180)),y);116 plot(U(:,floor(1+pi*40/(Dx*180)))/(2*Up*sin(40*pi/180)),y);117 plot(U(:,floor(1+pi*60/(Dx*180)))/(2*Up*sin(60*pi/180)),y);118 plot(U(:,floor(1+pi*80/(Dx*180)))/(2*Up*sin(80*pi/180)),y);119 plot(U(:,floor(1+pi*90/(Dx*180)))/(2*Up*sin(90*pi/180)),y);120 plot(U(:,floor(1+pi*100/(Dx*180)))/(2*Up*sin(100*pi/180)),y);121 plot(U(:,n)/(2*Up*sin(x(n))),y);122 axis([0,1.01,0,3])123 t=legend('$\varphi=20^\circ$','$\varphi=40^\circ$','$\varphi=60^\circ

$','$\varphi=80^\circ$','$\varphi=90^\circ$','$\varphi=100^\circ$','$\varphi=105.48^\circ$','Location','northwest');

124 set(t,'interpreter','latex')125

126 %Coeficiente de presión dependiente de x127 Cpx=zeros(Nx,1);128 Cpxteor=zeros(Nx,1);129 for i=1:Nx130 Cpx(i,1)=2*p(i);131 Cpxteor(i,1)=1−4*sin(x(i))^2;132 end133 figure134 plot(x,Cpx,'b')135 hold on136 plot(x,Cpxteor,'r')137 title('Coeficiente de presión en función de la posición')138 xlabel('$\alpha$','interpreter','latex')139 ylabel('$Cp(x)','interpreter','latex')140 grid on141 t=legend('Coeficiente de presión obtenido','Coeficiente de presión teórico'

,'location','southeast');142 set(t,'interpreter','latex')143 end

112 ETSII-UPM


C.3. Estenosis

C.3.1. Capa límite con condiciones de contorno de Poiseuille

1 %−−−−−−Programa cálculo de ley de d por ajuste por mínimos cuadrados−−−−−−%2 %d=exp(A)/Re^(1/B)3 %log(d)=A−1/B*log(Re)4 R1=3.5; %en mm5 x=20:10:2000; %Vector de números de Reynolds6 Re=[100 500 1000];7 d=[0.18 0.1 0.05]; %Valores de d estimados por iteraciones8 %Construcción del sistema M*incógnitas=b9 b=zeros(3,1);

10 b(:,1)=log(d);11 M=zeros(3,2);12 M(:,1)=1;13 M(:,2)=−log(Re);14 inc=M\b;15 A=inc(1,1);16 B=1/(inc(2,1));17 y=exp(A)./x.^(1/B); %d (es adimensional)18 fprintf('La ley para d es d=A/Re^(1/B) con B=%d\n',B)19 ylam=5./sqrt(x); %Adimensionalizado20 figure21 p1=loglog(x,y);22 hold on23 loglog(x(9),y(9),'ro') %Re=10024 hold on25 loglog(x(49),y(49),'ro') %Re=50026 hold on27 loglog(x(99),y(99),'ro') %Re=100028 hold on29 loglog(x(199),y(199),'go') %Re=200030 hold on31 p2=loglog(x,ylam);32 grid on33 xlabel('Re')34 ylabel('d')35 title('Ley \delta=f(Re) para flujo laminar')36 t=legend([p1,p2],'Flujo sanguineo','Placa plana');37 set(t,'interpreter','latex')

A modo de ejemplo se presenta el programa para resolver el caso M3. Para el resto decasos sería igual cambiando lo datos de entrada pertinentes

1 %−−−−−−−−−−Programa para la resolución del caso M3 (estenosis)−−−−−−−−−−%2 Reyn=[100 500 1000];3 tauups=zeros(length(Reyn),1);4 tauthroat=zeros(length(Reyn),1);



5 %Bucle que calcula para todos los Re6 for m=1:length(Reyn)7 %Datos8 R1=3.5*10^−3;9 ro=1060;

10 mu=3.71*10^−3;11 Re=Reyn(m);12 A=0.75;13 L=4*R1;14 %Datos derivados de los anteriores15 U1=mu*Re/(ro*R1);16 R2=R1*sqrt(1−A);17 if Re==10018 d=0.18;19 elseif Re==50020 d=0.1;21 elseif Re==100022 d=0.05;23 else24 d=0.0399;25 end26 U2=U1*R1^2/(2*(R2−d*R2)^2);27 long=sqrt((L/2)^2+(R1−R2)^2);28 %Adimensionalización de las variables29 ymax=10;30 U1a=ymax/sqrt(Re)*(2−ymax/sqrt(Re));31 U2a=U2/U1;32 %Contrucción de los vectores33 Dx=0.001;34 x=0:Dx:(2*long/R1);35 Dy=0.01;36 y=0:Dy:ymax;37 Nx=length(x);38 %Construcción del vector de velocidades y presiones39 Up=zeros(Nx,1);40 dp=zeros(Nx,1);41 for i=1:floor(Nx/2)42 Up(i)=(U2a−U1a)/(long/R1)*x(i)+U1a;43 %Up significa la de Poiseuille44 end45 k=i;46 for i=k:Nx47 Up(i)=−(U2a−U1a)/(long/R1)*(x(i)−x(k))+U2a;48 end49 Up=Up*L/(2*long);50 dU=diff(Up)/Dx;51 for i=1:(Nx−1)52 dp(i)=−Up(i)*dU(i);53 end54 dp(Nx)=dp(Nx−1);55 dp(floor(Nx/2)−1)=(dp(floor(Nx/2)−2)+dp(floor(Nx/2)))/2;56 %Resolución del problema57 N=length(y);58 U=zeros(N,Nx);

114 ETSII-UPM


59 V=zeros(N,Nx);60 uold=ones(N,1);61 for i=1:N62 uold(i)=y(i)/sqrt(Re)*(2−y(i)/sqrt(Re));63 end64 vold=zeros(N,1);65 U(:,1)=uold;66 V(:,1)=vold;67 A=zeros(N,N);68 B=zeros(N,1);69 K=zeros(N,N);70 L=zeros(N,1);71

72 %Bucle temporal;73 for n=2:Nx74 for i=2:N−175 A(i,i−1)=−vold(i)/Dy−1/Dy^2;76 A(i,i)=vold(i)/Dy+uold(i)/Dx+2/Dy^2;77 A(i,i+1)=−1/Dy^2;78 B(i,1)=−dp(n)+(uold(i))^2/Dx;79 end80 A(1,:)=0; A(1,1)=1;81 A(N,:)=0; A(N,N)=1;82 B(1,1)=0; B(N,1)=Up(n);83 Usol=A\B;84 Usol(1)=0;85 U(:,n)=Usol;86 for i=2:N87 K(i,i−1)=−1/Dy;88 K(i,i)=1/Dy;89 L(i,1)=−(U(i,n)−uold(i,1))/Dx;90 end91 K(1,:)=0; K(1,1)=1;92 L(1)=0;93 Vsol=K\L;94 Vsol(1)=0;95 V(:,n)=Vsol;96 uold=Usol;97 vold=Vsol;98 %Dejo de resolver cuando se desprende la capa límite99 if U(10,n)<0

100 break101 end102 end103 if U(10,n)<0104 fprintf('Se desprende en n=%d de Nx=%d\n',n,Nx)105 fprintf('Para Re=%d \n',Re)106 else107 n=Nx;108 end109 %Representacion de U110 figure111 plot(x(1)+U(:,1),y);112 title('U');



113 xlabel('Velocidad horizontal según la longitud adimensionalizada');114 ylabel('Altura');115 hold on116 grid on117 for i=1:n118 if mod(i,50)==0119 plot(x(i)+U(:,i),y)120 end121 end122 plot(x(n)+U(:,n),y)123 %Tau con dimensiones justo antes de la garganta (pico) (Upstream) y en la124 %gargante (throat)125 t=zeros(Nx,1);126 for i=1:Nx127 t(i)=mu*U1/R1*sqrt(Re)*(U(2,i)−U(1,i))/Dy*10;128 end129 tau_ups=max(t(Nx/4:Nx));130 tau_throat=t(floor(Nx/2));131 fprintf('Tau upstream vale %d dinas/cm^2\n',tau_ups)132 tauups(m)=tau_ups;133 tauthroat(m)=tau_throat;134 %Cálculo del caudal135 if ymax>(R2*sqrt(Re)/R1)136 K=0;137 [~,position]=min(abs(y−R2*sqrt(Re)/R1));138 for i=1:(position−1)139 K=K+Dy/2*(−R1/sqrt(Re)*(U(i,floor(Nx/2))*y(i)+U(i+1,floor(Nx/2))*y(

i+1))+R2*(U(i,floor(Nx/2))+U(i+1,floor(Nx/2))));140 end141 K=K*2*pi*U1*R1/sqrt(Re);142 Q2=K;143 else144 K=0;145 yq=0:Dy:(R2*sqrt(Re)/R1);146 U2q=ones(length(yq),1)*U2a;147 for i=1:N148 U2q(i,1)=U(i,floor(Nx/2));149 end150 for i=1:(length(yq)−1)151 K=K+Dy/2*(−R1/sqrt(Re)*(U2q(i,1)*yq(i)+U2q(i+1,1)*yq(i+1))+R2*(U2q(

i,1)+U2q(i+1,1)));152 end153 K=K*2*pi*U1*R1/sqrt(Re);154 Q2=K;155 end156 Q1=pi*R1^2*U1/2;157 fprintf('Caudal en 1 %d\n',Q1);158 fprintf('Caudal en 2 %d\n',Q2);159 error=abs(Q1−Q2)/Q1*100;160 fprintf('Error cometido en caudales %d\n',error);161 %Cálculo del espesor de desplazamiento162 K=0;163 Utr=U(:,floor(Nx/2));164 [~,position]=min(abs(Utr−0.999*U2a*4*R1/(2*long)));

116 ETSII-UPM


165 %Aquí he metido L=4R1166 for i=1:(position−1)167 K=K+Dy/2*(1−Utr(i)/(U2a*4*R1/(2*long))+1−Utr(i+1)/(U2a*4*R1/(2*long)));168 end169 delta1=K*R1/sqrt(Re);170 fprintf('Espesor de desplazamiento calculado %d y supuesto %d\n',delta1,d*

R2)171 end172 figure173 plot(Reyn,tauups,'r')174 hold on175 plot(Reyn,tauthroat,'b')176 hold on177 plot(Reyn,taupoi,'g')178 hold on179 plot(100,36.7,'v')180 hold on181 plot(500,316.9,'v')182 hold on183 plot(1000,788.1,'v')184 grid on185 legend('Upstream of throat','In throat','Poiseuille in throat')186 title('Variación de tau en función del Reynolds para M3')187 xlabel('Re')188 ylabel('tau en dinas/cm^2')189 save('tau_ups_M3','tauups')

C.3.2. Capa límite con condiciones de contorno de ujo potencial

1 //−−−−−−−−−Programa resolución flujo potencial en estenosis−−−−−−−−−//2 //Datos3 //Para los distintos casos hay que cambiar L, A y el nombre del archivo

de presiones4 real R1=0.0035*1000, A=0.44;5 real L=4*R1;6 real U1=0.5;7 real R2=R1*sqrt(1−A);8 real Dx=0.01;9

10 //Construcción del mallado11 border C1 (t=0,L/2) x=t;y=0;;12 border C2 (t=L/2,L) x=t;y=(R1−R2)*2/L*t+R2−R1;;13 border C3 (t=L,3*L/2) x=t;y=−(R1−R2)*2/L*t−3*(R2−R1);;14 border C4 (t=3*L/2,2*L) x=t;y=0;;15 border C5 (t=0,R1) x=2*L;y=t;;16 border C6 (t=2*L,0) x=t;y=R1;;17 border C7 (t=R1,0) x=0;y=t;;18 mesh Th=buildmesh(C1(10*L)+C2(20*L)+C3(20*L)+C4(10*L)+C5(10*R1)+C6(4*L*10)+

C7(10*R1));19 plot(Th,wait=0,cmm="Mallado");



20

21 //Resolución del problema22 fespace Vh2(Th,P2);23 Vh2 phi,g;24 fespace Vh1(Th,P1);25 Vh1 vx,vy,p;26 //Resolución de la ecuación lapl(phi)=027 problem laplace(phi,g)=int2d(Th)(dx(phi)*dx(g))+int2d(Th)(dy(phi)*dy(g))28 +on(C1,C2,C3,C4,phi=0)29 +on(C5,C7,phi=U1*y)30 +on(C6,phi=U1*R1);31

32 laplace;33 //Representación de la solución34 plot(phi,wait=0,nbiso=20,cmm=" Funcines de corriente ");35 vx=dy(phi);36 vy=−dx(phi);37 plot([vx,vy],wait=0,cmm=" Velocidades en x y en y");38 p =0.5*(1−vx^2 − vy^2);39 plot(p,ColorScheme=1,fill=1,nbiso=20,wait=0,cmm=" Campo de presiones ");40 //Lectura de la presión a una distancia h de la pared del vaso41 real alpha=atan(R1−R2,L/2);42 real h=0.12*R2;43 ofstream pr("presionesM1.txt");44 for (real i=0;i<=L/2;i=i+Dx*cos(alpha))45 pr<<p(L/2+i−h*sin(alpha),(R1−R2)*2/L*(i+L/2)+R2−R1+h*cos(alpha))<<endl;46 for (real i=0;i<=L/2;i=i+Dx*cos(alpha))47 pr<<p(i+L+h*sin(alpha),−(R1−R2)*2/L*(i+L)−3*(R2−R1)+h*cos(alpha))<<endl

;

1 %−−Programa resolución capa límite con cc de flujo potencial(estenosis)−−%2 %Todo está en uds del SI3 Dxff=0.01;4 R1=3.5*10^−3;5 p=load('presionesM1.txt');6 %long es la longitud total de la estenosis7 long=length(p)*Dxff*0.001;8 Nx=length(p)−1;9 Dx=long/(R1*Nx);

10 dp=diff(p)/Dx;11 x=0:Dx:(long/R1−Dx);12 Dy=0.01;13 ymax=10;14 y=0:Dy:ymax;15

16 figure17 plot(x,dp,'+−')18 grid on19

20 N=length(y);21 U=zeros(N,Nx);22 V=zeros(N,Nx);

118 ETSII-UPM


23 uold=ones(N,1);24 uold=0.5*uold;25 vold=zeros(N,1);26 U(:,1)=uold;27 V(:,1)=vold;28 A=zeros(N,N);29 B=zeros(N,1);30 K=zeros(N,N);31 L=zeros(N,1);32 %Bucle temporal;33 for n=2:Nx34 for i=2:N−135 A(i,i−1)=−vold(i)/Dy−1/Dy^2;36 A(i,i)=vold(i)/Dy+uold(i)/Dx+2/Dy^2;37 A(i,i+1)=−1/Dy^2;38 B(i,1)=−dp(n)+(uold(i))^2/Dx;39 end40 A(1,:)=0; A(1,1)=1;41 A(N,:)=0; A(N,N−1)=−1/Dy; A(N,N)=1/Dy;42 B(1,1)=0; B(N,1)=0;43 Usol=A\B;44 Usol(1)=0;45 U(:,n)=Usol;46 for i=2:N47 K(i,i−1)=−1/Dy;48 K(i,i)=1/Dy;49 L(i,1)=−(U(i,n)−uold(i,1))/Dx;50 end51 K(1,:)=0; K(1,1)=1;52 L(1)=0;53 Vsol=K\L;54 Vsol(1)=0;55 V(:,n)=Vsol;56 uold=Usol;57 vold=Vsol;58 %Dejo de resolver cuando se desprende la capa límite59 if U(10,n)<060 break61 end62 end63 if U(10,n)<064 fprintf('Se desprende en n=%d de Nx=%d\n',n,Nx)65 fprintf('Que es la longitud real en mm %d\n',x(n)*R1*1000)66 else67 n=Nx;68 end69 %Representacion de U70 figure71 plot(x(1)+U(:,1),y);72 title('U');73 xlabel('Velocidad horizontal según la longitud adimensionalizada');74 ylabel('Altura');75 hold on76 grid on



77 for i=1:n78 if mod(i,50)==079 plot(x(i)+U(:,i),y)80 end81 end82 plot(x(n)+U(:,n),y)

120 ETSII-UPM

Índice de guras

1. Flujo alrededor de un cuerpo . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2. Succión de la capa límite [12] . . . . . . . . . . . . . . . . . . . . . . . . . . 7

3. Perles de velocidades en distintos puntos de la placa plana . . . . . . . . . . 8

4. Geometría real y aproximada de la estenosis arterial . . . . . . . . . . . . . . 8

5. Resultados obtenidos: esfuerzo cortante máximo en la estenosis . . . . . . . . 9

1.1. Fuerzas de fricción y de presión sobre un perl aerodinámico . . . . . . . . . 4

1.2. Flujo alrededor de un cuerpo . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.1. Líneas de potencial constante y líneas de corriente ortogonales . . . . . . . . 12

2.2. Flujo ideal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.3. Flujo real . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3.1. Efecto del gradiente de presión adverso . . . . . . . . . . . . . . . . . . . . . 17

3.2. Capa límite adherida y desprendida . . . . . . . . . . . . . . . . . . . . . . . 18

3.3. Cuña de ángulo πβ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.4. Curvas de Falker-Skan. Fuente [12] . . . . . . . . . . . . . . . . . . . . . . . 21

3.5. Condición de contorno sobre f ′′ . . . . . . . . . . . . . . . . . . . . . . . . . 23

3.6. Curva obtenida, igual a las de Falker-Skan . . . . . . . . . . . . . . . . . . . 24

3.7. Perl de velocidades obtenido: perl de Blasius . . . . . . . . . . . . . . . . 25

4.1. Ejemplo de mallado en FreeFem++ . . . . . . . . . . . . . . . . . . . . . . 30

121

ÍNDICE DE FIGURAS

4.2. Mallado para los elementos n− 1, n y n+ 1 . . . . . . . . . . . . . . . . . . 33

5.1. Perl de entrada a la placa plana . . . . . . . . . . . . . . . . . . . . . . . . 40

5.2. Condiciones de contorno para ujo potencial en la placa plana . . . . . . . . 41

5.3. Perles de velocidades a lo largo de la placa . . . . . . . . . . . . . . . . . . 43

5.4. Comparación con las solución pseudoanalítica . . . . . . . . . . . . . . . . . 43

5.5. Coeciente de fricción en ejes logarítmicos . . . . . . . . . . . . . . . . . . . 44

5.6. Error cometido en Cf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

5.7. Resultados obtenidos Cf en la placa plana . . . . . . . . . . . . . . . . . . . 47

5.8. Succión en la capa límite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

5.9. Curvas para succión homogénea. Fuente [12] . . . . . . . . . . . . . . . . . . 48

5.10. Funciones de corriente en la placa plana para el caso de succión homogénea . 52

5.11. Comparación de los perles de velocidades. . . . . . . . . . . . . . . . . . . . 53

5.12. Distribución de velocidades en la placa placa con succión homogénea . . . . 54

5.13. Distribución de velocidades en la placa placa con succión homogénea, sin suc-ción y perl teórico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

6.1. Gradiente de presión alrededor de un cilindro . . . . . . . . . . . . . . . . . 56

6.2. Inuencia de la rugosidad en el carácter de la capa límite . . . . . . . . . . . 57

6.3. Mallado sobre el cilindro y funciones de corriente . . . . . . . . . . . . . . . 59

6.4. Valores de Cp en función de h . . . . . . . . . . . . . . . . . . . . . . . . . . 60

6.5. Gradiente de presiones a lo largo del cilindro . . . . . . . . . . . . . . . . . . 61

6.6. Perles de velocidades en el cilindro . . . . . . . . . . . . . . . . . . . . . . . 62

6.7. Ampliación del último perl de velocidades en el punto de desprendimiento. . 62

6.8. Comparación de las curvas obtenidas. Fuente [12] . . . . . . . . . . . . . . . 63

6.9. Coeciente de presiones en función de la posición . . . . . . . . . . . . . . . 64

6.10. Funciones de corriente con succión homogénea en el cilindro . . . . . . . . . 66

122 ETSII-UPM


6.11. Comparación de los perles de velocidades . . . . . . . . . . . . . . . . . . . 67

7.1. Angiografía de arteria carótida. Estrechamiento del 80% [6] . . . . . . . . . 70

7.2. Geometría de la estenosis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

7.3. Perles de velocidades en la estenosis . . . . . . . . . . . . . . . . . . . . . . 73

7.4. Espesor de desplazamiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

7.5. Variación de d con respecto a la variación de δ en el caso de una placa plana 78

7.6. Velocidad en los nodos superiores . . . . . . . . . . . . . . . . . . . . . . . . 80

7.7. Comparación de τmax calculado con los experimentales . . . . . . . . . . . . 81

7.8. Geometría del problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

7.9. Puntos de desprendimiento de la capa límite . . . . . . . . . . . . . . . . . . 87

7.10. Ampliaciones de los puntos de desprendimiento . . . . . . . . . . . . . . . . 87

7.11. Geometría de la estenosis más realista . . . . . . . . . . . . . . . . . . . . . . 89

A.1. EDP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

A.2. Diagrama de Gantt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97


ÍNDICE DE FIGURAS

124 ETSII-UPM

Índice de cuadros

6.1. Variación del ángulo de desprendimiento según la velocidad de succión . . . 66

7.1. Parámetros geométricos para el modelo . . . . . . . . . . . . . . . . . . . . . 72

7.2. Esfuerzo cortante máximo experimental (dina cm−2) . . . . . . . . . . . . . . 74

7.3. Parámetro de corrección d para los distintos Re . . . . . . . . . . . . . . . . 77

7.4. Comparación con los resultados experimentales y con los del artículo Reese(dina cm−2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

7.5. Caudales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

7.6. Espesor de desplazamiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

7.7. Espesor de la capa límite en la garganta . . . . . . . . . . . . . . . . . . . . 85

A.1. Presupuesto del proyecto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

125

ÍNDICE DE CUADROS

126 ETSII-UPM

Nomenclatura

δ Espesor de la capa límite

δ∗ Espesor de desplazamiento

µ Viscosidad dinámica del uido

ν Viscosidad cinemática del uido

φ Función potencial de la velocidad (Capítulo II)

ρ Densidad del uido

τ Esfuerzo cortante

ϕ Función de corriente

CD Coeciente de arrastre

Cf Coeciente de fricción

CL Coeciente de sustentación

Cp Coeciente de presión

D Longitud característica

Ff Fuerza de fricción

Fp Fuerza de presión

g Función base (Capítulo IV)

Nx Número de nodos en el eje X

Q Caudal

Re Número de Reynolds

u(x0, y) Perl de velocidades a la entrada de la geometría

U∞ Velocidad de la corriente en el innito

127

ÍNDICE DE CUADROS

Up Velocidad de la corriente de ujo potencial

v Velocidad

v0 Velocidad de succión

128 ETSII-UPM

Bibliografía

[1] Freefem++ website. http://www.freefem.org/.

[2] Matlab website. http://es.mathworks.com/products/matlab/index.html, 2016.

[3] Asai Asaithambi. A nite-dierence method for the falker-skan equation. 1998.

[4] L. H. Back & D. W. Crawford. Wall shear stress estimates in coronary artery constric-tions. 1992.

[5] Henry Rouviere & Andres Delmas. Anatomía Humana: descriptiva, topográca y fun-cional, pages 214217.

[6] David Sai Wah Ho et al. Multivessel spasm during coronary and peripheral angiography.2001.

[7] Richard T. Schoephoerster et al. Eects of local geometry ad uid dynamics on regionalplatelet deposition on articial surfaces. 1993.

[8] Rui Zhao et al. Micro-ow visualization of red blood cell-enhanced platelet concentrational sudden expansion. 2008.

[9] V. J. Modi & B. R. Seymour H.Huang. Fluid mechanics of stenoses arteries. 1995.

[10] David S. Thomson Jason M. Reese. Shear stress in arterial stenoses: a momentumintegral model. 1998.

[11] Antonio Crespo Martínez. Mecánica de Fluidos, pages 430436, 444446. 2006.

[12] Herrmann Schlichting. Teoría de la Capa Límite, pages 162177, 197, 378387. 1972.

[13] Theodoros G. Papaioannou & Chistodoulos Stefanadis. Vascular wall shear stress: Basicprinciples and methods. 2005.

[14] Frank M. White. Fluid Mechanics, chapter VII, VIII, pages 437563.

129

TRABAJO DE FIN DE GRADO - Archivo Digital UPM -...

Documents

Transcript of TRABAJO DE FIN DE GRADO - Archivo Digital UPM -...