MateWiki - Contribuciones del usuario [es]

Visor histórico del ferrocarril de Cuatro Vientos: análisis de la evolución del suelo y propuesta senda ciclable.

2024-05-22T12:32:17Z

Julia Cabeza:

{{ TrabajoSIG | Visor histórico del ferrocarril de Cuatro Vientos: análisis de la evolución del suelo y propuesta senda ciclable.
| Julia Cabeza Duque

Nuria Martínez Ballester

Sandra González de Mendoza Martínez|

[[:Categoría:SIGAIC_23/24|Curso 23/24]] }}

Este artículo recoge un trabajo de investigación, parte de la asignatura de [[:Categoría:Sistemas de Información Geográfica Aplicados a la Ingeniería Civil|Sistemas de Información Geográfica Aplicados a la Ingeniería Civil]].

En este artículo se va a realizar una visión a lo largo de los años del ferrocarril de Cuatro Vientos a Leganés. Históricamente era un ferrocarril de uso militar cuyo uso se ha visto opacado por la expansión de los núcleos urbanos próximos. El trabajo ha consistido en un análisis histórico del trazado del ferrocarril de Cuatro Vientos a Leganés. Para ello hemos empleado mapas e imágenes aéreas de distintos años (procedentes principalmente del IGN) y visualizado el crecimiento de los núcleos urbanos en torno a la línea.

Se ha comprobado la correcta georreferenciación de los mapas históricos y georreferenciado aquellos que obtuvimos en formato pdf o imágen. Para las diferentes series de años se ha obtenido el área de ferrocarril inicial que circula por área urbanizada y no urbanizada y se ha propuesto el trazado de un sendero ciclable con el que se aproveche el terreno inutilizado y mantenga el recuerdo de su pasado histórico.

== INTRODUCCIÓN ==

En los siguientes planos, obtenidos del servicio topográfico parcelario de 1941, podemos observar el carácter militar característico de Leganés en esta época. Se observan parcelas catalogadas como campos de instrucción y cuarteles, cuyos vestigios a día de hoy son prácticamente nulos: sólo quedan en los municipios próximos algunos centros deportivos militares. A pesar de que el carácter del suelo urbano ha cambiado drásticamente, su propiedad sigue perteneciendo al Ministerio de Defensa y como recuerdo de esa historia nos queda el trazado de la línea de ferrocarril de cuatro vientos.

[[Archivo: PARCELARIO 1 1.png|400px|marco|centro|Polígono 01 Leganés. Fuente: Catastro Histórico]]

[[Archivo: PARCELARIO 2.png|400px||marco|centro|Polígono 04 Leganés. Fuente: Catastro Histórico]]

[[Archivo: PARCELARIO 3.png|400px||marco|centro|Polígono 21 Leganés. Fuente: Catastro Histórico]]

Marcada en los anteriores planos como un área amarilla, se encuentra el área que ocupaba el ferrocarril mencionado durante su época auge. El trazado de vía ancha entre Cuatro Vientos y Leganés se abrió al servicio el 11 de enero de 1926, y servía principalmente como transporte de artillería. Este tramo fue construido con el ancho ibérico lo que permitía su conexión con el resto de la red ferroviaria. En 1943 se realiza una prolongación del trazado hasta Campamento.

El ferrocarril entra en desuso en 2002 debido a la rápida urbanización de los alrededores que hace que su función quede completamente obsoleta.

== METODOLOGÍA ==

Primero se ha realizado un estudio de la zona, obteniendo la localización de la línea y las estaciones de ferrocarril dentro de los municipios así como la topografía del terreno gracias a los modelos digitales del terreno descargados del CNIG.

El análisis es principalmente histórico. Se ha procedido a obtener mapas históricos de la zona en los que se observen los núcleos urbanos próximos y la evolución de la propia línea de ferrocarril. Se han empleado mapas, planos e imágenes aéreas de los siguientes años:

- 1875: previo a su construcción. Primera edición del MTN 50 del IGN mediante un servicio web conforme al estándar WMTS con dirección URL. Fuente: IGN

- 1937: previo a la prolongación a campamento. Hoja MTN 50 0559 Madrid de 1937 y Hoja MTN 50 0582 de Getafe. Fuente: IGN.

- 1941: parcelario. Parcelas 01, 04 y 21 de Leganés. Fuente: Catastro Histórico.

- 1956: Hojas MTN 50: 0582 de Getafe y 0559 de Madrid. Georreferenciadas. Fuente: IGN.

- 1975: con la estación en funcionamiento y recorrido final. Primera edición del MTN 25 del IGN mediante un servicio web conforme al estándar WMTS con dirección URL. Fuente: IGN

- 1998: se empieza a degradar su funcionamiento. PNOA huso 30 hoja 0559 Madrid y PNOA Huso 30 hoja 0582 Getafe del 1998. Fuente: IGN

- 2024: Imágenes aéreas, PNOA de máxima actualidad del IGN.

En algunos casos como el parcelario de 1941, el archivo fue obtenido en versión pdf por lo que hubo que georreferenciarlo para poder trabajar adecuadamente. Sirviendo como referencia el PNOA de máxima actualidad se procede con la herramienta Georeferencer de Ráster a su georeferenciación. Se emplearon 3 puntos y una georreferenciación lineal gracias a la buena calidad del archivo de origen.

A través de los años se ha clasificado el área de suelo que ocupaba inicialmente el ferrocarril como terreno urbanizado o no urbanizado con archivos vectoriales (obteniendo el área de cada tramo con la calculadora de campos de la tabla de atributos), observando la rapidez con la que los ayuntamientos avanzaban sobre el terreno militar, sobre todo en los últimos años.

Además, se ha realizado el posible trazado de una senda ciclable marcando los tramos peor conservados y los que necesitan menor rehabilitación. Se ha obtenido un perfil del trazado para verificar la viabilidad del sendero con el Plugin Profile Tool de QGIS.

== RESULTADOS ==

'''''LOCALIZACIÓN.'''''

Para la localización de la zona de trabajo obtuvimos las siluetas de los municipios de Madrid. Vemos que la línea discurría por Madrid, Alcorcón y Leganés.

[[Archivo:LOCALIZACION.png|500px|sinmarco|centro|Evolución de los principales núcleos urbanos]]

En concreto, para el estudio de la evolución del área urbanizada nos hemos enfocado en el tramo de Leganés, ya que estudios previos nos facilitaban esta información y planos como el de 1941 sólo señalaban este área:

[[Archivo:TRAMO ESTUDIADO.png|500px|sinmarco|centro|Evolución de los principales núcleos urbanos]]

Se obtuvieron las cotas de la zona directamente con el ráster inicial (MDT 25 Hojas 0559 y 0582). Se puede observar que el desnivel en toda la extensión es bastante bajo (cotas máximas de 700 m y mínimas de 670).

[[Archivo:ELEVACIONES.png|400px|sinmarco|centro|Evolución de los principales núcleos urbanos]]

'''''EVOLUCIÓN HISTÓRICA DE LOS NÚCLEOS URBANOS.'''''

En el siguiente mapa podemos observar la expansión de los núcleos urbanos próximos a la línea (los núcleos urbanos próximos más importantes son Leganés al sudeste, Carabanchel al norte y Alcorcón al sudoeste). Esta expansión es la causante de la pérdida por parte del ejército no de la propiedad de sus terrenos si no del uso sobre el área, es decir el suelo sigue siendo suyo pero se utiliza para vivienda principalmente. Además, observamos el trazado de dos líneas, una previa a la extensión de la línea de 1943, y la final que sí que incluye el tramo hasta Campamento.

[[Archivo:CONJUNTO NUCLEOS URBANOS.png|750px|marco|centro|Evolución de los principales núcleos urbanos]]

'''''EVOLUCIÓN SUELO URBANO.'''''

En las siguientes imágenes se observa la serie de 1875-2024 del suelo no urbanizado y el urbanizado. En concreto se trata del tramo de Leganés (señalado en el apartado de Introducción) que en principio abarcaba un área no urbanizada de 180,593.62 m2 (1937), actualmente el área en contacto con zona no urbanizada se limita a unos 35,656.02 m2:

[[Archivo:CONJUNTO AREA URBANIZADA 3.png|750px|marco|centro|Evolución área urbanizada en torno al ferrocarril]]

'''''SENDA CICLABLE.'''''

Una propuesta de actuación sobre el trazado sería la creación de un sendero ciclable de 8km de extensión recorriendo el antiguo trazado de la línea de ferrocarril. Se han señalado las zonas peor conservadas que requieren de una actuación mayor. Además, se propone la instalación de parques saludables (con ejercicios aptos para la tercera edad) al inicio y final de tramo. Además del carril bici, se adaptará un sendero adyacente peatonal separado con vegetación para la seguridad de los peatones.

[[Archivo:CARRIL BICI COMPLETO.png|marco|centro|Senda ciclabe propuesta]]

[[Archivo:PERFIL CARRIL.png|marco|centro|Perfil del carril bici propuesto]]

== CONCLUSIONES ==
Como se ha podido analizar, el suelo urbano correspondiente a los distintos ayuntamientos se ha
expandido notablemente hacia las zonas de propiedad militar. Esta expansión no se ha realizado de
forma legal. El suelo sigue siendo propiedad del Ministerio de Defensa, pero este no hace uso del
mismo.

Esta situación genera irregularidades legales: si ocurre algún accidente en una de las parcelas
edificadas debido a una construcción negligente o por anomalías en el terreno, el responsable último
sería la Administración (Ministerio de Defensa).

La solución legal a este problema sería la cesión por parte de los ayuntamientos de los municipios de
una parcela edificable de valor igual al área que se ha ocupado. La propuesta de la senda ciclable
permitiría que cada ayuntamiento hiciese uso del tramo que le correspondiese, otorgando al
Ministerio una parcela equivalente.

== ANEXO ==

En los siguientes apartados recogemos todos los planos con detalle que hemos generado en QGIS
para el mejor entendimiento del análisis realizado.

'''''EVOLUCIÓN HISTÓRICA.'''''

-Mapa de los núcleos urbanos en 1875

[[Archivo:MAPA NUCLEOS 1875.png|750px|sinmarco|centro|Mapa de los núcleos urbanos en 1875]]

-Trazado futura línea de ferrocarril en mapa 1875

[[Archivo:FUTURA LINEA FERROCARRIL.png|750px|sinmarco|centro|Trazado futura línea de ferrocarril en mapa 1875]]

-Mapa de los núcleos urbanos en 1937

[[Archivo:MAPA NUCLEOS 1937.png|750px|sinmarco|centro|Mapa de los núcleos urbanos en 1937]]

-Mapa de los núcleos urbanos en 1956

[[Archivo:MAPA NUCLEOS 1956.png|750px|sinmarco|centro|Mapa de los núcleos urbanos en 1956]]

-Mapa de los núcleos urbanos en 1975

[[Archivo:MAPA NUCLEOS 1975.png|750px|sinmarco|centro|Mapa de los núcleos urbanos en 1975]]

-Mapa de los núcleos urbanos en 1998

[[Archivo:NUCLEOS URBANOS 1998.png|750px|sinmarco|centro|Mapa de los núcleos urbanos en 1998]]

-Mapa de los núcleos urbanos ACTUALES

[[Archivo:MAPA NUCLEOS ACTUAL .png|750px|sinmarco|centro|Mapa de los núcleos urbanos ACTUALES]]

'''''EVOLUCIÓN DEL SUELO URBANO.'''''

-División área urbanizada frente área no urbanizada de 1875

[[Archivo:AREA URBANIZADA 1875.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada de 1875]]

-División área urbanizada frente área no urbanizada de 1937

[[Archivo:AREA URBANIZADA 1937.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada de 1937]]

-División área urbanizada frente área no urbanizada de 1956

[[Archivo:AREA URBANIZADA 1956.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada 1956]]

-División área urbanizada frente área no urbanizada de 1975

[[Archivo:AREA URBANIZADA 1975.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada de 1975]]

-División área urbanizada frente área no urbanizada de 1998

[[Archivo:AREA URBANIZADA 1998.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada de 1998]]

-División área urbanizada frente área no urbanizada ACTUAL

[[Archivo:AREA URBANIZADA ACTUAL.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada ACTUAL]]

'''''PROPUESTA SENDERO.'''''

-Detalle del carril bici trazado:

[[Archivo:DETALLE CARRIL 1.png|750px|sinmarco|centro]]

[[Archivo:DETALLE CARRIL 2.png|750px|sinmarco|centro]]

[[Archivo:DETALLE CARRIL 3.png|750px|sinmarco|centro]]

[[Archivo:DETALLE CARRIL 4.png|750px|sinmarco|centro]]

[[Categoría:Sistemas de Información Geográfica Aplicados a la Ingeniería Civil]]
[[Categoría:SIGAIC_23/24]]

Visor histórico del ferrocarril de Cuatro Vientos: análisis de la evolución del suelo y propuesta senda ciclable.

2024-05-22T12:30:24Z

Julia Cabeza:

{{ TrabajoSIG | Visor histórico del ferrocarril de Cuatro Vientos: análisis de la evolución del suelo y propuesta senda ciclable.
| Julia Cabeza Duque

Nuria Martínez Ballester

Sandra González de Mendoza Martínez|

[[:Categoría:SIGAIC_23/24|Curso 23/24]] }}

Este artículo recoge un trabajo de investigación, parte de la asignatura de [[:Categoría:Sistemas de Información Geográfica Aplicados a la Ingeniería Civil|Sistemas de Información Geográfica Aplicados a la Ingeniería Civil]].

En este artículo se va a realizar una visión a lo largo de los años del ferrocarril de Cuatro Vientos a Leganés. Históricamente era un ferrocarril de uso militar cuyo uso se ha visto opacado por la expansión de los núcleos urbanos próximos. El trabajo ha consistido en un análisis histórico del trazado del ferrocarril de Cuatro Vientos a Leganés. Para ello hemos empleado mapas e imágenes aéreas de distintos años (procedentes principalmente del IGN) y visualizado el crecimiento de los núcleos urbanos en torno a la línea.

Se ha comprobado la correcta georreferenciación de los mapas históricos y georreferenciado aquellos que obtuvimos en formato pdf o imágen. Para las diferentes series de años se ha obtenido el área de ferrocarril inicial que circula por área urbanizada y no urbanizada y se ha propuesto el trazado de un sendero ciclable con el que se aproveche el terreno inutilizado y mantenga el recuerdo de su pasado histórico.

== INTRODUCCIÓN ==

En los siguientes planos, obtenidos del servicio topográfico parcelario de 1941, podemos observar el carácter militar característico de Leganés en esta época. Se observan parcelas catalogadas como campos de instrucción y cuarteles, cuyos vestigios a día de hoy son prácticamente nulos: sólo quedan en los municipios próximos algunos centros deportivos militares. A pesar de que el carácter del suelo urbano ha cambiado drásticamente, su propiedad sigue perteneciendo al Ministerio de Defensa y como recuerdo de esa historia nos queda el trazado de la línea de ferrocarril de cuatro vientos.

[[Archivo: PARCELARIO 1 1.png|400px|marco|centro|Polígono 01 Leganés. Fuente: Catastro Histórico]]

[[Archivo: PARCELARIO 2.png|400px||marco|centro|Polígono 04 Leganés. Fuente: Catastro Histórico]]

[[Archivo: PARCELARIO 3.png|400px||marco|centro|Polígono 21 Leganés. Fuente: Catastro Histórico]]

Marcada en los anteriores planos como un área amarilla, se encuentra el área que ocupaba el ferrocarril mencionado durante su época auge. El trazado de vía ancha entre Cuatro Vientos y Leganés se abrió al servicio el 11 de enero de 1926, y servía principalmente como transporte de artillería. Este tramo fue construido con el ancho ibérico lo que permitía su conexión con el resto de la red ferroviaria. En 1943 se realiza una prolongación del trazado hasta Campamento.

El ferrocarril entra en desuso en 2002 debido a la rápida urbanización de los alrededores que hace que su función quede completamente obsoleta.

== METODOLOGÍA ==

Primero se ha realizado un estudio de la zona, obteniendo la localización de la línea y las estaciones de ferrocarril dentro de los municipios así como la topografía del terreno gracias a los modelos digitales del terreno descargados del CNIG.

El análisis es principalmente histórico. Se ha procedido a obtener mapas históricos de la zona en los que se observen los núcleos urbanos próximos y la evolución de la propia línea de ferrocarril. Se han empleado mapas, planos e imágenes aéreas de los siguientes años:

- 1875: previo a su construcción. Primera edición del MTN 50 del IGN mediante un servicio web conforme al estándar WMTS con dirección URL. Fuente: IGN

- 1937: previo a la prolongación a campamento. Hoja MTN 50 0559 Madrid de 1937 y Hoja MTN 50 0582 de Getafe. Fuente: IGN.

- 1941: parcelario. Parcelas 01, 04 y 21 de Leganés. Fuente: Catastro Histórico.

- 1956: Hojas MTN 50: 0582 de Getafe y 0559 de Madrid. Georreferenciadas. Fuente: IGN.

- 1975: con la estación en funcionamiento y recorrido final. Primera edición del MTN 25 del IGN mediante un servicio web conforme al estándar WMTS con dirección URL. Fuente: IGN

- 1998: se empieza a degradar su funcionamiento. PNOA huso 30 hoja 0559 Madrid y PNOA Huso 30 hoja 0582 Getafe del 1998. Fuente: IGN

- 2024: Imágenes aéreas, PNOA de máxima actualidad del IGN.

En algunos casos como el parcelario de 1941, el archivo fue obtenido en versión pdf por lo que hubo que georreferenciarlo para poder trabajar adecuadamente. Sirviendo como referencia el PNOA de máxima actualidad se procede con la herramienta Georeferencer de Ráster a su georeferenciación. Se emplearon 3 puntos y una georreferenciación lineal gracias a la buena calidad del archivo de origen.

A través de los años se ha clasificado el área de suelo que ocupaba inicialmente el ferrocarril como terreno urbanizado o no urbanizado con archivos vectoriales (obteniendo el área de cada tramo con la calculadora de campos de la tabla de atributos), observando la rapidez con la que los ayuntamientos avanzaban sobre el terreno militar, sobre todo en los últimos años.

Además, se ha realizado el posible trazado de una senda ciclable marcando los tramos peor conservados y los que necesitan menor rehabilitación. Se ha obtenido un perfil del trazado para verificar la viabilidad del sendero con el Plugin Profile Tool de QGIS.

== RESULTADOS ==

'''''LOCALIZACIÓN.'''''

Para la localización de la zona de trabajo obtuvimos las siluetas de los municipios de Madrid. Vemos que la línea discurría por Madrid, Alcorcón y Leganés.

[[Archivo:LOCALIZACION.png|500px|sinmarco|centro|Evolución de los principales núcleos urbanos]]

En concreto, para el estudio de la evolución del área urbanizada nos hemos enfocado en el tramo de Leganés, ya que estudios previos nos facilitaban esta información y planos como el de 1941 sólo señalaban este área:

[[Archivo:TRAMO ESTUDIADO.png|500px|sinmarco|centro|Evolución de los principales núcleos urbanos]]

Se obtuvieron las cotas de la zona directamente con el ráster inicial (MDT 25 Hojas 0559 y 0582). Se puede observar que el desnivel en toda la extensión es bastante bajo (cotas máximas de 700 m y mínimas de 670).

[[Archivo:ELEVACIONES.png|400px|sinmarco|centro|Evolución de los principales núcleos urbanos]]

'''''EVOLUCIÓN HISTÓRICA DE LOS NÚCLEOS URBANOS.'''''

En el siguiente mapa podemos observar la expansión de los núcleos urbanos próximos a la línea (los núcleos urbanos próximos más importantes son Leganés al sudeste, Carabanchel al norte y Alcorcón al sudoeste). Esta expansión es la causante de la pérdida por parte del ejército no de la propiedad de sus terrenos si no del uso sobre el área, es decir el suelo sigue siendo suyo pero se utiliza para vivienda principalmente. Además, observamos el trazado de dos líneas, una previa a la extensión de la línea de 1943, y la final que sí que incluye el tramo hasta Campamento.

[[Archivo:CONJUNTO NUCLEOS URBANOS.png|750px|marco|centro|Evolución de los principales núcleos urbanos]]

'''''EVOLUCIÓN SUELO URBANO.'''''

En las siguientes imágenes se observa la serie de 1875-2024 del suelo no urbanizado y el urbanizado. En concreto se trata del tramo de Leganés (señalado en el apartado de Introducción) que en principio abarcaba un área no urbanizada de 180,593.62 m2 (1937), actualmente el área en contacto con zona no urbanizada se limita a unos 35,656.02 m2:

[[Archivo:CONJUNTO AREA URBANIZADA 3.png|750px|marco|centro|Evolución área urbanizada en torno al ferrocarril]]

'''''SENDA CICLABLE.'''''

Una propuesta de actuación sobre el trazado sería la creación de un sendero ciclable de 8km de extensión recorriendo el antiguo trazado de la línea de ferrocarril. Se han señalado las zonas peor conservadas que requieren de una actuación mayor. Además, se propone la instalación de parques saludables (con ejercicios aptos para la tercera edad) al inicio y final de tramo. Además del carril bici, se adaptará un sendero adyacente peatonal separado con vegetación para la seguridad de los peatones.

[[Archivo:CARRIL BICI COMPLETO.png|marco|centro|Senda ciclabe propuesta]]

[[Archivo:PERFIL CARRIL.png|marco|centro|Perfil del carril bici propuesto]]

== CONCLUSIONES ==
Como se ha podido analizar, el suelo urbano correspondiente a los distintos ayuntamientos se ha
expandido notablemente hacia las zonas de propiedad militar. Esta expansión no se ha realizado de
forma legal. El suelo sigue siendo propiedad del Ministerio de Defensa, pero este no hace uso del
mismo.

Esta situación genera irregularidades legales: si ocurre algún accidente en una de las parcelas
edificadas debido a una construcción negligente o por anomalías en el terreno, el responsable último
sería la Administración (Ministerio de Defensa).

La solución legal a este problema sería la cesión por parte de los ayuntamientos de los municipios de
una parcela edificable de valor igual al área que se ha ocupado. La propuesta de la senda ciclable
permitiría que cada ayuntamiento hiciese uso del tramo que le correspondiese, otorgando al
Ministerio una parcela equivalente.

== ANEXO ==

En los siguientes apartados recogemos todos los planos con detalle que hemos generado en QGIS
para el mejor entendimiento del análisis realizado.

'''''EVOLUCIÓN HISTÓRICA.'''''

-Mapa de los núcleos urbanos en 1875

[[Archivo:MAPA NUCLEOS 1875.png|750px|sinmarco|centro|Mapa de los núcleos urbanos en 1875]]

-Trazado futura línea de ferrocarril en mapa 1875

[[Archivo:FUTURA LINEA FERROCARRIL.png|750px|sinmarco|centro|Trazado futura línea de ferrocarril en mapa 1875]]

-Mapa de los núcleos urbanos en 1937

[[Archivo:MAPA NUCLEOS 1937.png|750px|sinmarco|centro|Mapa de los núcleos urbanos en 1937]]

-Mapa de los núcleos urbanos en 1956

[[Archivo:MAPA NUCLEOS 1956.png|750px|sinmarco|centro|Mapa de los núcleos urbanos en 1956]]

-Mapa de los núcleos urbanos en 1975

[[Archivo:MAPA NUCLEOS 1975.png|750px|sinmarco|centro|Mapa de los núcleos urbanos en 1975]]

-Mapa de los núcleos urbanos en 1998

[[Archivo:NUCLEOS URBANOS 1998.png|750px|sinmarco|centro|Mapa de los núcleos urbanos en 1998]]

-Mapa de los núcleos urbanos ACTUALES

[[Archivo:MAPA NUCLEOS ACTUAL .png|750px|sinmarco|centro|Mapa de los núcleos urbanos ACTUALES]]

'''''EVOLUCIÓN DEL SUELO URBANO.'''''

-División área urbanizada frente área no urbanizada de 1875

[[Archivo:AREA URBANIZADA 1875.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada de 1875]]

-División área urbanizada frente área no urbanizada de 1937

[[Archivo:AREA URBANIZADA 1937.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada de 1937]]

-División área urbanizada frente área no urbanizada de 1956

[[Archivo:AREA URBANIZADA 1956.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada 1956]]

-División área urbanizada frente área no urbanizada de 1975

[[Archivo:AREA URBANIZADA 1975.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada de 1975]]

-División área urbanizada frente área no urbanizada de 1998

[[Archivo:AREA URBANIZADA 1998.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada de 1998]]

-División área urbanizada frente área no urbanizada ACTUAL

[[Archivo:AREA URBANIZADA ACTUAL.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada ACTUAL]]

'''''PROPUESTA SENDERO.'''''

-Detalle del carril bici trazado:

[[Archivo:DETALLE CARRIL 1.png|750px|sinmarco|centro]]

[[Archivo:DETALLE CARRIL 2.png|750px|sinmarco|centro]]

[[Archivo:DETALLE CARRIL 3.png|750px|sinmarco|centro]]

[[Archivo:DETALLE CARRIL 4.png|750px|sinmarco|centro]]

<pre><nowiki>[[Categoría:SIGAIC_23/24]]</nowiki></pre>

Visor histórico del ferrocarril de Cuatro Vientos: análisis de la evolución del suelo y propuesta senda ciclable.

2024-05-19T13:18:54Z

Julia Cabeza:

{{ TrabajoSIG | Visor histórico del ferrocarril de Cuatro Vientos: análisis de la evolución del suelo y propuesta senda ciclable| Julia Cabeza Duque, Nuria Martínez Ballester, Sandra González de Mendoza Martínez.| [[:Categoría:SIGAIC_23/24|Curso 23/24]] }}

En este artículo se va a realizar una visión a lo largo de los años del ferrocarril de Cuatro Vientos a Leganés. Históricamente era un ferrocarril de uso militar cuyo uso se ha visto opacado por la expansión de los núcleos urbanos próximos. El trabajo ha consistido en un análisis histórico del trazado del ferrocarril de Cuatro Vientos a Leganés. Para ello hemos empleado mapas e imágenes aéreas de distintos años (procedentes principalmente del IGN) y visualizado el crecimiento de los núcleos urbanos en torno a la línea.

Se ha comprobado la correcta georreferenciación de los mapas históricos y georreferenciado aquellos que obtuvimos en formato pdf o imágen. Para las diferentes series de años se ha obtenido el área de ferrocarril inicial que circula por área urbanizada y no urbanizada y se ha propuesto el trazado de un sendero ciclable con el que se aproveche el terreno inutilizado y mantenga el recuerdo de su pasado histórico.

== INTRODUCCIÓN ==

En los siguientes planos, obtenidos del servicio topográfico parcelario de 1941, podemos observar el carácter militar característico de Leganés en esta época. Se observan parcelas catalogadas como campos de instrucción y cuarteles, cuyos vestigios a día de hoy son prácticamente nulos: sólo quedan en los municipios próximos algunos centros deportivos militares. A pesar de que el carácter del suelo urbano ha cambiado drásticamente, su propiedad sigue perteneciendo al Ministerio de Defensa y como recuerdo de esa historia nos queda el trazado de la línea de ferrocarril de cuatro vientos.

[[Archivo: PARCELARIO 1 1.png|400px|marco|centro|Polígono 01 Leganés. Fuente: Catastro Histórico]]

[[Archivo: PARCELARIO 2.png|400px||marco|centro|Polígono 04 Leganés. Fuente: Catastro Histórico]]

[[Archivo: PARCELARIO 3.png|400px||marco|centro|Polígono 21 Leganés. Fuente: Catastro Histórico]]

Marcada en los anteriores planos como un área amarilla, se encuentra el área que ocupaba el ferrocarril mencionado durante su época auge. El trazado de vía ancha entre Cuatro Vientos y Leganés se abrió al servicio el 11 de enero de 1926, y servía principalmente como transporte de artillería. Este tramo fue construido con el ancho ibérico lo que permitía su conexión con el resto de la red ferroviaria. En 1943 se realiza una prolongación del trazado hasta Campamento.

El ferrocarril entra en desuso en 2002 debido a la rápida urbanización de los alrededores que hace que su función quede completamente obsoleta.

== METODOLOGÍA ==

Primero se ha realizado un estudio de la zona, obteniendo la localización de la línea y las estaciones de ferrocarril dentro de los municipios así como la topografía del terreno gracias a los modelos digitales del terreno descargados del CNIG.

El análisis es principalmente histórico. Se ha procedido a obtener mapas históricos de la zona en los que se observen los núcleos urbanos próximos y la evolución de la propia línea de ferrocarril. Se han empleado mapas, planos e imágenes aéreas de los siguientes años:

- 1875: previo a su construcción. Primera edición del MTN 50 del IGN mediante un servicio web conforme al estándar WMTS con dirección URL. Fuente: IGN

- 1937: previo a la prolongación a campamento. Hoja MTN 50 0559 Madrid de 1937 y Hoja MTN 50 0582 de Getafe. Fuente: IGN.

- 1941: parcelario. Parcelas 01, 04 y 21 de Leganés. Fuente: Catastro Histórico.

- 1956: Hojas MTN 50: 0582 de Getafe y 0559 de Madrid. Georreferenciadas. Fuente: IGN.

- 1975: con la estación en funcionamiento y recorrido final. Primera edición del MTN 25 del IGN mediante un servicio web conforme al estándar WMTS con dirección URL. Fuente: IGN

- 1998: se empieza a degradar su funcionamiento. PNOA huso 30 hoja 0559 Madrid y PNOA Huso 30 hoja 0582 Getafe del 1998. Fuente: IGN

- 2024: Imágenes aéreas, PNOA de máxima actualidad del IGN.

En algunos casos como el parcelario de 1941, el archivo fue obtenido en versión pdf por lo que hubo que georreferenciarlo para poder trabajar adecuadamente. Sirviendo como referencia el PNOA de máxima actualidad se procede con la herramienta Georeferencer de Ráster a su georeferenciación. Se emplearon 3 puntos y una georreferenciación lineal gracias a la buena calidad del archivo de origen.

A través de los años se ha clasificado el área de suelo que ocupaba inicialmente el ferrocarril como terreno urbanizado o no urbanizado con archivos vectoriales (obteniendo el área de cada tramo con la calculadora de campos de la tabla de atributos), observando la rapidez con la que los ayuntamientos avanzaban sobre el terreno militar, sobre todo en los últimos años.

Además, se ha realizado el posible trazado de una senda ciclable marcando los tramos peor conservados y los que necesitan menor rehabilitación. Se ha obtenido un perfil del trazado para verificar la viabilidad del sendero con el Plugin Profile Tool de QGIS.

== RESULTADOS ==

'''''LOCALIZACIÓN.'''''

Para la localización de la zona de trabajo obtuvimos las siluetas de los municipios de Madrid. Vemos que la línea discurría por Madrid, Alcorcón y Leganés.

[[Archivo:LOCALIZACION.png|500px|sinmarco|centro|Evolución de los principales núcleos urbanos]]

En concreto, para el estudio de la evolución del área urbanizada nos hemos enfocado en el tramo de Leganés, ya que estudios previos nos facilitaban esta información y planos como el de 1941 sólo señalaban este área:

[[Archivo:TRAMO ESTUDIADO.png|500px|sinmarco|centro|Evolución de los principales núcleos urbanos]]

Se obtuvieron las cotas de la zona directamente con el ráster inicial (MDT 25 Hojas 0559 y 0582). Se puede observar que el desnivel en toda la extensión es bastante bajo (cotas máximas de 700 m y mínimas de 670).

[[Archivo:ELEVACIONES.png|400px|sinmarco|centro|Evolución de los principales núcleos urbanos]]

'''''EVOLUCIÓN HISTÓRICA DE LOS NÚCLEOS URBANOS.'''''

En el siguiente mapa podemos observar la expansión de los núcleos urbanos próximos a la línea (los núcleos urbanos próximos más importantes son Leganés al sudeste, Carabanchel al norte y Alcorcón al sudoeste). Esta expansión es la causante de la pérdida por parte del ejército no de la propiedad de sus terrenos si no del uso sobre el área, es decir el suelo sigue siendo suyo pero se utiliza para vivienda principalmente. Además, observamos el trazado de dos líneas, una previa a la extensión de la línea de 1943, y la final que sí que incluye el tramo hasta Campamento.

[[Archivo:CONJUNTO NUCLEOS URBANOS.png|750px|marco|centro|Evolución de los principales núcleos urbanos]]

'''''EVOLUCIÓN SUELO URBANO.'''''

En las siguientes imágenes se observa la serie de 1875-2024 del suelo no urbanizado y el urbanizado. En concreto se trata del tramo de Leganés (señalado en el apartado de Introducción) que en principio abarcaba un área no urbanizada de 180,593.62 m2 (1937), actualmente el área en contacto con zona no urbanizada se limita a unos 35,656.02 m2:

[[Archivo:CONJUNTO AREA URBANIZADA 3.png|750px|marco|centro|Evolución área urbanizada en torno al ferrocarril]]

'''''SENDA CICLABLE.'''''

Una propuesta de actuación sobre el trazado sería la creación de un sendero ciclable de 8km de extensión recorriendo el antiguo trazado de la línea de ferrocarril. Se han señalado las zonas peor conservadas que requieren de una actuación mayor. Además, se propone la instalación de parques saludables (con ejercicios aptos para la tercera edad) al inicio y final de tramo. Además del carril bici, se adaptará un sendero adyacente peatonal separado con vegetación para la seguridad de los peatones.

[[Archivo:CARRIL BICI COMPLETO.png|marco|centro|Senda ciclabe propuesta]]

[[Archivo:PERFIL CARRIL.png|marco|centro|Perfil del carril bici propuesto]]

== CONCLUSIONES ==
Como se ha podido analizar, el suelo urbano correspondiente a los distintos ayuntamientos se ha
expandido notablemente hacia las zonas de propiedad militar. Esta expansión no se ha realizado de
forma legal. El suelo sigue siendo propiedad del Ministerio de Defensa, pero este no hace uso del
mismo.

Esta situación genera irregularidades legales: si ocurre algún accidente en una de las parcelas
edificadas debido a una construcción negligente o por anomalías en el terreno, el responsable último
sería la Administración (Ministerio de Defensa).

La solución legal a este problema sería la cesión por parte de los ayuntamientos de los municipios de
una parcela edificable de valor igual al área que se ha ocupado. La propuesta de la senda ciclable
permitiría que cada ayuntamiento hiciese uso del tramo que le correspondiese, otorgando al
Ministerio una parcela equivalente.

== ANEXO ==

En los siguientes apartados recogemos todos los planos con detalle que hemos generado en QGIS
para el mejor entendimiento del análisis realizado.

'''''EVOLUCIÓN HISTÓRICA.'''''

-Mapa de los núcleos urbanos en 1875

[[Archivo:MAPA NUCLEOS 1875.png|750px|sinmarco|centro|Mapa de los núcleos urbanos en 1875]]

-Trazado futura línea de ferrocarril en mapa 1875

[[Archivo:FUTURA LINEA FERROCARRIL.png|750px|sinmarco|centro|Trazado futura línea de ferrocarril en mapa 1875]]

-Mapa de los núcleos urbanos en 1937

[[Archivo:MAPA NUCLEOS 1937.png|750px|sinmarco|centro|Mapa de los núcleos urbanos en 1937]]

-Mapa de los núcleos urbanos en 1956

[[Archivo:MAPA NUCLEOS 1956.png|750px|sinmarco|centro|Mapa de los núcleos urbanos en 1956]]

-Mapa de los núcleos urbanos en 1975

[[Archivo:MAPA NUCLEOS 1975.png|750px|sinmarco|centro|Mapa de los núcleos urbanos en 1975]]

-Mapa de los núcleos urbanos en 1998

[[Archivo:NUCLEOS URBANOS 1998.png|750px|sinmarco|centro|Mapa de los núcleos urbanos en 1998]]

-Mapa de los núcleos urbanos ACTUALES

[[Archivo:MAPA NUCLEOS ACTUAL .png|750px|sinmarco|centro|Mapa de los núcleos urbanos ACTUALES]]

'''''EVOLUCIÓN DEL SUELO URBANO.'''''

-División área urbanizada frente área no urbanizada de 1875

[[Archivo:AREA URBANIZADA 1875.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada de 1875]]

-División área urbanizada frente área no urbanizada de 1937

[[Archivo:AREA URBANIZADA 1937.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada de 1937]]

-División área urbanizada frente área no urbanizada de 1956

[[Archivo:AREA URBANIZADA 1956.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada 1956]]

-División área urbanizada frente área no urbanizada de 1975

[[Archivo:AREA URBANIZADA 1975.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada de 1975]]

-División área urbanizada frente área no urbanizada de 1998

[[Archivo:AREA URBANIZADA 1998.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada de 1998]]

-División área urbanizada frente área no urbanizada ACTUAL

[[Archivo:AREA URBANIZADA ACTUAL.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada ACTUAL]]

'''''PROPUESTA SENDERO.'''''

-Detalle del carril bici trazado:

[[Archivo:DETALLE CARRIL 1.png|750px|sinmarco|centro]]

[[Archivo:DETALLE CARRIL 2.png|750px|sinmarco|centro]]

[[Archivo:DETALLE CARRIL 3.png|750px|sinmarco|centro]]

[[Archivo:DETALLE CARRIL 4.png|750px|sinmarco|centro]]

Archivo:DETALLE CARRIL 4.png

2024-05-19T13:17:49Z

Julia Cabeza:

Archivo:DETALLE CARRIL 3.png

2024-05-19T13:17:11Z

Julia Cabeza:

Archivo:DETALLE CARRIL 2.png

2024-05-19T13:14:43Z

Julia Cabeza:

Archivo:DETALLE CARRIL 1.png

2024-05-19T13:14:16Z

Julia Cabeza:

Visor histórico del ferrocarril de Cuatro Vientos: análisis de la evolución del suelo y propuesta senda ciclable.

2024-05-19T12:01:27Z

Julia Cabeza:

Archivo:MAPA NUCLEOS 1975.png

2024-05-19T11:59:28Z

Julia Cabeza:

Archivo:MAPA NUCLEOS 1937.png

2024-05-19T11:58:25Z

Julia Cabeza:

Archivo:MAPA NUCLEOS ACTUAL .png

2024-05-19T11:57:32Z

Julia Cabeza:

Archivo:MAPA NUCLEOS 1956.png

2024-05-19T11:56:48Z

Julia Cabeza:

Archivo:FUTURA LINEA FERROCARRIL.png

2024-05-19T11:52:59Z

Julia Cabeza:

Archivo:PERFIL CARRIL.png

2024-05-19T11:50:41Z

Julia Cabeza:

Archivo:CARRIL BICI COMPLETO.png

2024-05-19T11:49:34Z

Julia Cabeza:

Archivo:AREA URBANIZADA 1998.png

2024-05-19T11:46:44Z

Julia Cabeza:

Archivo:NUCLEOS URBANOS 1998.png

2024-05-19T11:44:34Z

Julia Cabeza:

Archivo:AREA URBANIZADA 1956.png

2024-05-19T11:41:22Z

Julia Cabeza:

Archivo:PARCELARIO 3.png

2024-05-19T11:38:38Z

Julia Cabeza:

Archivo:PARCELARIO 2.png

2024-05-19T11:38:12Z

Julia Cabeza:

Archivo:PARCELARIO 1 1.png

2024-05-19T11:37:14Z

Julia Cabeza:

Archivo:PARCELARIO 1.png

2024-05-19T11:32:14Z

Julia Cabeza:

Visor histórico del ferrocarril de Cuatro Vientos: análisis de la evolución del suelo y propuesta senda ciclable.

2024-05-18T18:00:57Z

Julia Cabeza:

{{ TrabajoSIG | Visor histórico del ferrocarril de Cuatro Vientos: análisis de la evolución del suelo y propuesta senda ciclable| Julia Cabeza Duque, Nuria Martínez Ballester, Sandra González de Mendoza Martínez.| [[:Categoría:SIGAIC_23/24|Curso 23/24]] }}

En este artículo se va a realizar una visión a lo largo de los años del ferrocarril de Cuatro Vientos a Leganés. Históricamente era un ferrocarril de uso militar cuyo uso se ha visto opacado por la expansión de los núcleos urbanos próximos. El trabajo ha consistido en un análisis histórico del trazado del ferrocarril de Cuatro Vientos a Leganés. Para ello hemos empleado mapas e imágenes aéreas de distintos años (procedentes principalmente del IGN) y visualizado el crecimiento de los núcleos urbanos en torno a la línea.

Se ha comprobado la correcta georreferenciación de los mapas históricos y georreferenciado aquellos que obtuvimos en formato pdf o imágen. Para las diferentes series de años se ha obtenido el área de ferrocarril inicial que circula por área urbanizada y no urbanizada y se ha propuesto el trazado de un sendero ciclable con el que se aproveche el terreno inutilizado y mantenga el recuerdo de su pasado histórico.

== INTRODUCCIÓN ==

En los siguientes planos, obtenidos del servicio topográfico parcelario de 1941, podemos observar el carácter militar característico de Leganés en esta época. Se observan parcelas catalogadas como campos de instrucción y cuarteles, cuyos vestigios a día de hoy son prácticamente nulos: sólo quedan en los municipios próximos algunos centros deportivos militares. A pesar de que el carácter del suelo urbano ha cambiado drásticamente, su propiedad sigue perteneciendo al Ministerio de Defensa y como recuerdo de esa historia nos queda el trazado de la línea de ferrocarril de cuatro vientos.

[[Archivo:SIG Cabeza 1.png|marco|centro|Polígono 04 Leganés. Fuente: Catastro Histórico]]

Marcada en los anteriores planos como un área amarilla, se encuentra el área que ocupaba el ferrocarril mencionado durante su época auge. El trazado de vía ancha entre Cuatro Vientos y Leganés se abrió al servicio el 11 de enero de 1926, y servía principalmente como transporte de artillería. Este tramo fue construido con el ancho ibérico lo que permitía su conexión con el resto de la red ferroviaria. En 1943 se realiza una prolongación del trazado hasta Campamento.

El ferrocarril entra en desuso en 2002 debido a la rápida urbanización de los alrededores que hace que su función quede completamente obsoleta.

== METODOLOGÍA ==

Primero se ha realizado un estudio de la zona, obteniendo la localización de la línea y las estaciones de ferrocarril dentro de los municipios así como la topografía del terreno gracias a los modelos digitales del terreno descargados del CNIG.

El análisis es principalmente histórico. Se ha procedido a obtener mapas históricos de la zona en los que se observen los núcleos urbanos próximos y la evolución de la propia línea de ferrocarril. Se han empleado mapas, planos e imágenes aéreas de los siguientes años:

- 1875: previo a su construcción. Primera edición del MTN 50 del IGN mediante un servicio web conforme al estándar WMTS con dirección URL. Fuente: IGN

- 1937: previo a la prolongación a campamento. Hoja MTN 50 0559 Madrid de 1937 y Hoja MTN 50 0582 de Getafe. Fuente: IGN.

- 1941: parcelario. Parcelas 01, 04 y 21 de Leganés. Fuente: Catastro Histórico.

- 1956: Hojas MTN 50: 0582 de Getafe y 0559 de Madrid. Georreferenciadas. Fuente: IGN.

- 1975: con la estación en funcionamiento y recorrido final. Primera edición del MTN 25 del IGN mediante un servicio web conforme al estándar WMTS con dirección URL. Fuente: IGN

- 1998: se empieza a degradar su funcionamiento. PNOA huso 30 hoja 0559 Madrid y PNOA Huso 30 hoja 0582 Getafe del 1998. Fuente: IGN

- 2024: Imágenes aéreas, PNOA de máxima actualidad del IGN.

En algunos casos como el parcelario de 1941, el archivo fue obtenido en versión pdf por lo que hubo que georreferenciarlo para poder trabajar adecuadamente. Sirviendo como referencia el PNOA de máxima actualidad se procede con la herramienta Georeferencer de Ráster a su georeferenciación. Se emplearon 3 puntos y una georreferenciación lineal gracias a la buena calidad del archivo de origen.

A través de los años se ha clasificado el área de suelo que ocupaba inicialmente el ferrocarril como terreno urbanizado o no urbanizado con archivos vectoriales (obteniendo el área de cada tramo con la calculadora de campos de la tabla de atributos), observando la rapidez con la que los ayuntamientos avanzaban sobre el terreno militar, sobre todo en los últimos años.

Además, se ha realizado el posible trazado de una senda ciclable marcando los tramos peor conservados y los que necesitan menor rehabilitación. Se ha obtenido un perfil del trazado para verificar la viabilidad del sendero con el Plugin Profile Tool de QGIS.

== RESULTADOS ==

'''''LOCALIZACIÓN.'''''

Para la localización de la zona de trabajo obtuvimos las siluetas de los municipios de Madrid. Vemos que la línea discurría por Madrid, Alcorcón y Leganés.

[[Archivo:LOCALIZACION.png|500px|sinmarco|centro|Evolución de los principales núcleos urbanos]]

En concreto, para el estudio de la evolución del área urbanizada nos hemos enfocado en el tramo de Leganés, ya que estudios previos nos facilitaban esta información y planos como el de 1941 sólo señalaban este área:

[[Archivo:TRAMO ESTUDIADO.png|500px|sinmarco|centro|Evolución de los principales núcleos urbanos]]

Se obtuvieron las cotas de la zona directamente con el ráster inicial (MDT 25 Hojas 0559 y 0582). Se puede observar que el desnivel en toda la extensión es bastante bajo (cotas máximas de 700 m y mínimas de 670).

[[Archivo:ELEVACIONES.png|400px|sinmarco|centro|Evolución de los principales núcleos urbanos]]

'''''EVOLUCIÓN HISTÓRICA DE LOS NÚCLEOS URBANOS.'''''

En el siguiente mapa podemos observar la expansión de los núcleos urbanos próximos a la línea (los núcleos urbanos próximos más importantes son Leganés al sudeste, Carabanchel al norte y Alcorcón al sudoeste). Esta expansión es la causante de la pérdida por parte del ejército no de la propiedad de sus terrenos si no del uso sobre el área, es decir el suelo sigue siendo suyo pero se utiliza para vivienda principalmente. Además, observamos el trazado de dos líneas, una previa a la extensión de la línea de 1943, y la final que sí que incluye el tramo hasta Campamento.

[[Archivo:CONJUNTO NUCLEOS URBANOS.png|750px|marco|centro|Evolución de los principales núcleos urbanos]]

'''''EVOLUCIÓN SUELO URBANO.'''''

En las siguientes imágenes se observa la serie de 1875-2024 del suelo no urbanizado y el urbanizado. En concreto se trata del tramo de Leganés (señalado en el apartado de Introducción) que en principio abarcaba un área no urbanizada de 180,593.62 m2 (1937), actualmente el área en contacto con zona no urbanizada se limita a unos 35,656.02 m2:

[[Archivo:CONJUNTO AREA URBANIZADA 3.png|750px|marco|centro|Evolución área urbanizada en torno al ferrocarril]]

'''''SENDA CICLABLE.'''''

Una propuesta de actuación sobre el trazado sería la creación de un sendero ciclable de 8km de extensión recorriendo el antiguo trazado de la línea de ferrocarril. Se han señalado las zonas peor conservadas que requieren de una actuación mayor. Además, se propone la instalación de parques saludables (con ejercicios aptos para la tercera edad) al inicio y final de tramo. Además del carril bici, se adaptará un sendero adyacente peatonal separado con vegetación para la seguridad de los peatones.

== CONCLUSIONES ==
Como se ha podido analizar, el suelo urbano correspondiente a los distintos ayuntamientos se ha
expandido notablemente hacia las zonas de propiedad militar. Esta expansión no se ha realizado de
forma legal. El suelo sigue siendo propiedad del Ministerio de Defensa, pero este no hace uso del
mismo.

Esta situación genera irregularidades legales: si ocurre algún accidente en una de las parcelas
edificadas debido a una construcción negligente o por anomalías en el terreno, el responsable último
sería la Administración (Ministerio de Defensa).

La solución legal a este problema sería la cesión por parte de los ayuntamientos de los municipios de
una parcela edificable de valor igual al área que se ha ocupado. La propuesta de la senda ciclable
permitiría que cada ayuntamiento hiciese uso del tramo que le correspondiese, otorgando al
Ministerio una parcela equivalente.

== ANEXO ==

En los siguientes apartados recogemos todos los planos con detalle que hemos generado en QGIS
para el mejor entendimiento del análisis realizado.

'''''EVOLUCIÓN HISTÓRICA.'''''

-1875

[[Archivo:MAPA NUCLEOS 1875.png|750px|sinmarco|centro|Mapa de los núcleos urbanos en 1875]]

'''''EVOLUCIÓN DEL SUELO URBANO.'''''

-1875

[[Archivo:AREA URBANIZADA 1875.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada de 1875]]

-1937

[[Archivo:AREA URBANIZADA 1937.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada de 1937]]

-1975

[[Archivo:AREA URBANIZADA 1975.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada de 1975]]

-ACTUAL

[[Archivo:AREA URBANIZADA ACTUAL.png|750px|sinmarco|centro| División área urbanizada frente área no urbanizada ACTUAL]]

'''''PROPUESTA SENDERO.'''''

Archivo:ELEVACIONES.png

2024-05-18T17:58:01Z

Julia Cabeza:

Archivo:TRAMO ESTUDIADO.png

2024-05-18T17:57:28Z

Julia Cabeza:

Archivo:LOCALIZACION.png

2024-05-18T17:56:27Z

Julia Cabeza:

Archivo:AREA URBANIZADA ACTUAL.png

2024-05-18T17:54:00Z

Julia Cabeza:

Archivo:CONJUNTO AREA URBANIZADA 3.png

2024-05-18T17:50:49Z

Julia Cabeza:

Archivo:CONJUNTO NUCLEOS URBANOS.png

2024-05-18T17:48:54Z

Julia Cabeza:

Archivo:CONJUNTO AREA URBANIZADA 2.png

2024-05-18T17:44:22Z

Julia Cabeza:

Archivo:CONJUNTO AREA URBANIZADA.png

2024-05-18T17:36:01Z

Julia Cabeza:

Archivo:AREA URBANIZADA 1937.png

2024-05-18T17:34:11Z

Julia Cabeza:

Visor histórico del ferrocarril de Cuatro Vientos: análisis de la evolución del suelo y propuesta senda ciclable.

2024-05-18T16:10:22Z

Julia Cabeza:

{{ TrabajoSIG | Visor histórico del ferrocarril de Cuatro Vientos: análisis de la evolución del suelo y propuesta senda ciclable| Julia Cabeza Duque, Nuria Martínez Ballester, Sandra González de Mendoza Martínez.| [[:Categoría:SIGAIC_23/24|Curso 23/24]] }}

En este artículo se va a realizar una visión a lo largo de los años del ferrocarril de Cuatro Vientos a Leganés. Históricamente era un ferrocarril de uso militar cuyo uso se ha visto opacado por la expansión de los núcleos urbanos próximos. El trabajo ha consistido en un análisis histórico del trazado del ferrocarril de Cuatro Vientos a Leganés. Para ello hemos empleado mapas e imágenes aéreas de distintos años (procedentes principalmente del IGN) y visualizado el crecimiento de los núcleos urbanos en torno a la línea.

Se ha comprobado la correcta georreferenciación de los mapas históricos y georreferenciado aquellos que obtuvimos en formato pdf o imágen. Para las diferentes series de años se ha obtenido el área de ferrocarril inicial que circula por área urbanizada y no urbanizada y se ha propuesto el trazado de un sendero ciclable con el que se aproveche el terreno inutilizado y mantenga el recuerdo de su pasado histórico.

== INTRODUCCIÓN ==

En los siguientes planos, obtenidos del servicio topográfico parcelario de 1941, podemos observar el carácter militar característico de Leganés en esta época. Se observan parcelas catalogadas como campos de instrucción y cuarteles, cuyos vestigios a día de hoy son prácticamente nulos: sólo quedan en los municipios próximos algunos centros deportivos militares. A pesar de que el carácter del suelo urbano ha cambiado drásticamente, su propiedad sigue perteneciendo al Ministerio de Defensa y como recuerdo de esa historia nos queda el trazado de la línea de ferrocarril de cuatro vientos.
[[Archivo:SIG Cabeza 1.png|750px|marco|centro|Polígono 04 Leganés. Fuente: Catastro Histórico]]

Marcada en los anteriores planos como un área amarilla, se encuentra el área que ocupaba el ferrocarril mencionado durante su época auge. El trazado de vía ancha entre Cuatro Vientos y Leganés se abrió al servicio el 11 de enero de 1926, y servía principalmente como transporte de artillería. Este tramo fue construido con el ancho ibérico lo que permitía su conexión con el resto de la red ferroviaria. En 1943 se realiza una prolongación del trazado hasta Campamento.

El ferrocarril entra en desuso en 2002 debido a la rápida urbanización de los alrededores que hace que su función quede completamente obsoleta.

== METODOLOGÍA ==

Primero se ha realizado un estudio de la zona, obteniendo la localización de la línea y las estaciones de ferrocarril dentro de los municipios así como la topografía del terreno gracias a los modelos digitales del terreno descargados del CNIG.

El análisis es principalmente histórico. Se ha procedido a obtener mapas históricos de la zona en los que se observen los núcleos urbanos próximos y la evolución de la propia línea de ferrocarril. Se han empleado mapas, planos e imágenes aéreas de los siguientes años:

- 1875: previo a su construcción. Primera edición del MTN 50 del IGN mediante un servicio web conforme al estándar WMTS con dirección URL. Fuente: IGN
- 1937: previo a la prolongación a campamento. Hoja MTN 50 0559 Madrid de 1937 y Hoja MTN 50 0582 de Getafe. Fuente: IGN.
- 1941: parcelario. Parcelas 01, 04 y 21 de Leganés. Fuente: Catastro Histórico.
- 1956: Hojas MTN 50: 0582 de Getafe y 0559 de Madrid. Georreferenciadas. Fuente: IGN.
- 1975: con la estación en funcionamiento y recorrido final. Primera edición del MTN 25 del IGN mediante un servicio web conforme al estándar WMTS con dirección URL. Fuente: IGN
- 1998: se empieza a degradar su funcionamiento. PNOA huso 30 hoja 0559 Madrid y PNOA Huso 30 hoja 0582 Getafe del 1998. Fuente: IGN
- 2024: Imágenes aéreas, PNOA de máxima actualidad del IGN.

En algunos casos como el parcelario de 1941, el archivo fue obtenido en versión pdf por lo que hubo que georreferenciarlo para poder trabajar adecuadamente. Sirviendo como referencia el PNOA de máxima actualidad se procede con la herramienta Georeferencer de Ráster a su georeferenciación. Se emplearon 3 puntos y una georreferenciación lineal gracias a la buena calidad del archivo de origen.

A través de los años se ha clasificado el área de suelo que ocupaba inicialmente el ferrocarril como terreno urbanizado o no urbanizado con archivos vectoriales (obteniendo el área de cada tramo con la calculadora de campos de la tabla de atributos), observando la rapidez con la que los ayuntamientos avanzaban sobre el terreno militar, sobre todo en los últimos años.

Además, se ha realizado el posible trazado de una senda ciclable marcando los tramos peor conservados y los que necesitan menor rehabilitación. Se ha obtenido un perfil del trazado para verificar la viabilidad del sendero con el Plugin Profile Tool de QGIS.

== RESULTADOS ==

'''''LOCALIZACIÓN.'''''

Para la localización de la zona de trabajo obtuvimos las siluetas de los municipios de Madrid. Vemos que la línea discurría por Madrid, Alcorcón y Leganés.

En concreto, para el estudio de la evolución del área urbanizada nos hemos enfocado en el tramo de Leganés, ya que estudios previos nos facilitaban esta información y planos como el de 1941 sólo señalaban este área:

Se obtuvieron las cotas de la zona directamente con el ráster inicial (MDT 25 Hojas 0559 y 0582). Se puede observar que el desnivel en toda la extensión es bastante bajo (cotas máximas de 700 m y mínimas de 670).

'''''EVOLUCIÓN HISTÓRICA DE LOS NÚCLEOS URBANOS.'''''

En el siguiente mapa podemos observar la expansión de los núcleos urbanos próximos a la línea (los núcleos urbanos próximos más importantes son Leganés al sudeste, Carabanchel al norte y Alcorcón al sudoeste). Esta expansión es la causante de la pérdida por parte del ejército no de la propiedad de sus terrenos si no del uso sobre el área, es decir el suelo sigue siendo suyo pero se utiliza para vivienda principalmente. Además, observamos el trazado de dos líneas, una previa a la extensión de la línea de 1943, y la final que sí que incluye el tramo hasta Campamento.

'''''EVOLUCIÓN SUELO URBANO.'''''

En las siguientes imágenes se observa la serie de 1875-2024 del suelo no urbanizado y el urbanizado. En concreto se trata del tramo de Leganés (señalado en el apartado de Introducción) que en principio abarcaba un área no urbanizada de 180,593.62 m2 (1937), actualmente el área en contacto con zona no urbanizada se limita a unos 35,656.02 m2:

'''''SENDA CICLABLE.'''''

Una propuesta de actuación sobre el trazado sería la creación de un sendero ciclable de 8km de extensión recorriendo el antiguo trazado de la línea de ferrocarril. Se han señalado las zonas peor conservadas que requieren de una actuación mayor. Además, se propone la instalación de parques saludables (con ejercicios aptos para la tercera edad) al inicio y final de tramo. Además del carril bici, se adaptará un sendero adyacente peatonal separado con vegetación para la seguridad de los peatones.

== CONCLUSIONES ==
Como se ha podido analizar, el suelo urbano correspondiente a los distintos ayuntamientos se ha
expandido notablemente hacia las zonas de propiedad militar. Esta expansión no se ha realizado de
forma legal. El suelo sigue siendo propiedad del Ministerio de Defensa, pero este no hace uso del
mismo.

Esta situación genera irregularidades legales: si ocurre algún accidente en una de las parcelas
edificadas debido a una construcción negligente o por anomalías en el terreno, el responsable último
sería la Administración (Ministerio de Defensa).

La solución legal a este problema sería la cesión por parte de los ayuntamientos de los municipios de
una parcela edificable de valor igual al área que se ha ocupado. La propuesta de la senda ciclable
permitiría que cada ayuntamiento hiciese uso del tramo que le correspondiese, otorgando al
Ministerio una parcela equivalente.

== ANEXO ==

En los siguientes apartados recogemos todos los planos con detalle que hemos generado en QGIS
para el mejor entendimiento del análisis realizado.

'''''EVOLUCIÓN HISTÓRICA.'''''

-1875

[[Archivo:MAPA NUCLEOS 1875.png|750px|sinmarco|centro|Mapa de los núcleos urbanos en 1875]]

'''''EVOLUCIÓN DEL SUELO URBANO.'''''

'''''PROPUESTA SENDERO.'''''

Archivo:MAPA NUCLEOS 1875.png

2024-05-18T16:03:01Z

Julia Cabeza:

Archivo:AREA URBANIZADA 1875.png

2024-05-18T16:00:32Z

Julia Cabeza:

Archivo:AREA URBANIZADA 1975.png

2024-05-18T15:59:24Z

Julia Cabeza:

Visor histórico del ferrocarril de Cuatro Vientos: análisis de la evolución del suelo y propuesta senda ciclable.

2024-05-18T15:58:11Z

Julia Cabeza:

{{ TrabajoSIG | Visor histórico del ferrocarril de Cuatro Vientos: análisis de la evolución del suelo y propuesta senda ciclable| Julia Cabeza Duque, Nuria Martínez Ballester, Sandra González de Mendoza Martínez.| [[:Categoría:SIGAIC_23/24|Curso 23/24]] }}

En este artículo se va a realizar una visión a lo largo de los años del ferrocarril de Cuatro Vientos a Leganés. Históricamente era un ferrocarril de uso militar cuyo uso se ha visto opacado por la expansión de los núcleos urbanos próximos. El trabajo ha consistido en un análisis histórico del trazado del ferrocarril de Cuatro Vientos a Leganés. Para ello hemos empleado mapas e imágenes aéreas de distintos años (procedentes principalmente del IGN) y visualizado el crecimiento de los núcleos urbanos en torno a la línea.

Se ha comprobado la correcta georreferenciación de los mapas históricos y georreferenciado aquellos que obtuvimos en formato pdf o imágen. Para las diferentes series de años se ha obtenido el área de ferrocarril inicial que circula por área urbanizada y no urbanizada y se ha propuesto el trazado de un sendero ciclable con el que se aproveche el terreno inutilizado y mantenga el recuerdo de su pasado histórico.

== INTRODUCCIÓN ==

En los siguientes planos, obtenidos del servicio topográfico parcelario de 1941, podemos observar el carácter militar característico de Leganés en esta época. Se observan parcelas catalogadas como campos de instrucción y cuarteles, cuyos vestigios a día de hoy son prácticamente nulos: sólo quedan en los municipios próximos algunos centros deportivos militares. A pesar de que el carácter del suelo urbano ha cambiado drásticamente, su propiedad sigue perteneciendo al Ministerio de Defensa y como recuerdo de esa historia nos queda el trazado de la línea de ferrocarril de cuatro vientos.
[[Archivo:SIG Cabeza 1.png|750px|marco|centro|Polígono 04 Leganés. Fuente: Catastro Histórico]]

Marcada en los anteriores planos como un área amarilla, se encuentra el área que ocupaba el ferrocarril mencionado durante su época auge. El trazado de vía ancha entre Cuatro Vientos y Leganés se abrió al servicio el 11 de enero de 1926, y servía principalmente como transporte de artillería. Este tramo fue construido con el ancho ibérico lo que permitía su conexión con el resto de la red ferroviaria. En 1943 se realiza una prolongación del trazado hasta Campamento.

El ferrocarril entra en desuso en 2002 debido a la rápida urbanización de los alrededores que hace que su función quede completamente obsoleta.

== METODOLOGÍA ==

Primero se ha realizado un estudio de la zona, obteniendo la localización de la línea y las estaciones de ferrocarril dentro de los municipios así como la topografía del terreno gracias a los modelos digitales del terreno descargados del CNIG.

El análisis es principalmente histórico. Se ha procedido a obtener mapas históricos de la zona en los que se observen los núcleos urbanos próximos y la evolución de la propia línea de ferrocarril. Se han empleado mapas, planos e imágenes aéreas de los siguientes años:

- 1875: previo a su construcción. Primera edición del MTN 50 del IGN mediante un servicio web conforme al estándar WMTS con dirección URL. Fuente: IGN
- 1937: previo a la prolongación a campamento. Hoja MTN 50 0559 Madrid de 1937 y Hoja MTN 50 0582 de Getafe. Fuente: IGN.
- 1941: parcelario. Parcelas 01, 04 y 21 de Leganés. Fuente: Catastro Histórico.
- 1956: Hojas MTN 50: 0582 de Getafe y 0559 de Madrid. Georreferenciadas. Fuente: IGN.
- 1975: con la estación en funcionamiento y recorrido final. Primera edición del MTN 25 del IGN mediante un servicio web conforme al estándar WMTS con dirección URL. Fuente: IGN
- 1998: se empieza a degradar su funcionamiento. PNOA huso 30 hoja 0559 Madrid y PNOA Huso 30 hoja 0582 Getafe del 1998. Fuente: IGN
- 2024: Imágenes aéreas, PNOA de máxima actualidad del IGN.

En algunos casos como el parcelario de 1941, el archivo fue obtenido en versión pdf por lo que hubo que georreferenciarlo para poder trabajar adecuadamente. Sirviendo como referencia el PNOA de máxima actualidad se procede con la herramienta Georeferencer de Ráster a su georeferenciación. Se emplearon 3 puntos y una georreferenciación lineal gracias a la buena calidad del archivo de origen.

A través de los años se ha clasificado el área de suelo que ocupaba inicialmente el ferrocarril como terreno urbanizado o no urbanizado con archivos vectoriales (obteniendo el área de cada tramo con la calculadora de campos de la tabla de atributos), observando la rapidez con la que los ayuntamientos avanzaban sobre el terreno militar, sobre todo en los últimos años.

Además, se ha realizado el posible trazado de una senda ciclable marcando los tramos peor conservados y los que necesitan menor rehabilitación. Se ha obtenido un perfil del trazado para verificar la viabilidad del sendero con el Plugin Profile Tool de QGIS.

== RESULTADOS ==

'''''LOCALIZACIÓN.'''''

Para la localización de la zona de trabajo obtuvimos las siluetas de los municipios de Madrid. Vemos que la línea discurría por Madrid, Alcorcón y Leganés.

En concreto, para el estudio de la evolución del área urbanizada nos hemos enfocado en el tramo de Leganés, ya que estudios previos nos facilitaban esta información y planos como el de 1941 sólo señalaban este área:

Se obtuvieron las cotas de la zona directamente con el ráster inicial (MDT 25 Hojas 0559 y 0582). Se puede observar que el desnivel en toda la extensión es bastante bajo (cotas máximas de 700 m y mínimas de 670).

'''''EVOLUCIÓN HISTÓRICA DE LOS NÚCLEOS URBANOS.'''''

En el siguiente mapa podemos observar la expansión de los núcleos urbanos próximos a la línea (los núcleos urbanos próximos más importantes son Leganés al sudeste, Carabanchel al norte y Alcorcón al sudoeste). Esta expansión es la causante de la pérdida por parte del ejército no de la propiedad de sus terrenos si no del uso sobre el área, es decir el suelo sigue siendo suyo pero se utiliza para vivienda principalmente. Además, observamos el trazado de dos líneas, una previa a la extensión de la línea de 1943, y la final que sí que incluye el tramo hasta Campamento.

'''''EVOLUCIÓN SUELO URBANO.'''''

En las siguientes imágenes se observa la serie de 1875-2024 del suelo no urbanizado y el urbanizado. En concreto se trata del tramo de Leganés (señalado en el apartado de Introducción) que en principio abarcaba un área no urbanizada de 180,593.62 m2 (1937), actualmente el área en contacto con zona no urbanizada se limita a unos 35,656.02 m2:

'''''SENDA CICLABLE.'''''

Una propuesta de actuación sobre el trazado sería la creación de un sendero ciclable de 8km de extensión recorriendo el antiguo trazado de la línea de ferrocarril. Se han señalado las zonas peor conservadas que requieren de una actuación mayor. Además, se propone la instalación de parques saludables (con ejercicios aptos para la tercera edad) al inicio y final de tramo. Además del carril bici, se adaptará un sendero adyacente peatonal separado con vegetación para la seguridad de los peatones.

== CONCLUSIONES ==
Como se ha podido analizar, el suelo urbano correspondiente a los distintos ayuntamientos se ha
expandido notablemente hacia las zonas de propiedad militar. Esta expansión no se ha realizado de
forma legal. El suelo sigue siendo propiedad del Ministerio de Defensa, pero este no hace uso del
mismo.

Esta situación genera irregularidades legales: si ocurre algún accidente en una de las parcelas
edificadas debido a una construcción negligente o por anomalías en el terreno, el responsable último
sería la Administración (Ministerio de Defensa).

La solución legal a este problema sería la cesión por parte de los ayuntamientos de los municipios de
una parcela edificable de valor igual al área que se ha ocupado. La propuesta de la senda ciclable
permitiría que cada ayuntamiento hiciese uso del tramo que le correspondiese, otorgando al
Ministerio una parcela equivalente.

== ANEXO ==

En los siguientes apartados recogemos todos los planos con detalle que hemos generado en QGIS
para el mejor entendimiento del análisis realizado.

'''''EVOLUCIÓN HISTÓRICA.'''''

'''''EVOLUCIÓN DEL SUELO URBANO.'''''

'''''PROPUESTA SENDERO.'''''

Archivo:SIG Cabeza 1.png

2024-05-18T15:33:03Z

Julia Cabeza:

Visor histórico del ferrocarril de Cuatro Vientos: análisis de la evolución del suelo y propuesta senda ciclable.

2024-05-12T13:22:13Z

Julia Cabeza:

Visor histórico del ferrocarril de Cuatro Vientos: análisis de la evolución del suelo y propuesta senda ciclable.

2024-05-12T13:20:51Z

Julia Cabeza:

Visor histórico del ferrocarril de Cuatro Vientos: análisis de la evolución del suelo y propuesta senda ciclable.

2024-05-12T13:11:01Z

Julia Cabeza: Página creada con «{{ TrabajoSIG | Visor histórico del ferrocarril de Cuatro Vientos: análisis de la evolución del suelo y propuesta senda ciclable| Julia Cabeza Duque, Nuria Martínez B...»

Visualización de campos escalares y vectoriales en fluidos. (Grupo B6)

2021-12-14T07:41:43Z

Julia Cabeza: /* º Puntos frontera */

{{ TrabajoED | Visualización de campos escalares y vectoriales en fluidos. (Grupo B6) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC21/22|2021-22]] | María Magaldi Jurado, Nuria Martínez Ballester, Julia Cabeza Duque }}

En este trabajo vamos a visualizar campos escalares y vectoriales en fluidos. Más específicamente, el fluido a tratar es uno incompresible que se encuentra alrededor de un obstáculo de forma circular, formando así un anillo comprendido entre las circunferencias de radio 1 y 5, con el centro en el origen de coordenadas. Para facilitar el trabajo usaremos coordenadas polares (cilíndricas <math>(\rho,\theta)</math>). Para dibujar que el fluido va por el exterior de la superficie definida consideraremos que <math>(\rho,\theta) \in [0,5]*[0,2\pi]</math>.
Además, tendremos que la velocidad de las partículas del fluido viene definida por el gradiente de la función potencial <math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math>

== º Mallado del fluido ==
El siguiente mallado representa los puntos interiores de la región ocupada por un fluido. Se encuentra en el exterior del círculo unidad. Como ya se ha introducido, nuestro mallado del anillo está comprendido entre los radios 1 y 5 y centro el origen.

[[Archivo:1,mallad0.jpg|800px|thumb|center|Mallado del fluido]]

El CÓDIGO matlab empleado es:
{{matlab|codigo=
ro=linspace(1,5,50); %creación de los intervalos
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt); %creación matrices
Mx=Mro.*cos(Mtt); %parametrización de x
My=Mro.*sin(Mtt); %parametrización de y
Mz=zeros(size(Mx)); %parametrización de z (nula al estar en un plano)
subplot(1,1,1);
mesh(Mx,My,Mz) %mallado
axis([-5,5,-5,5]) %superficie definida en [-5,5]*[-5,5]
view(2)
}}

== º Función potencial <math>\varphi</math> y velocidad de las partículas del fluido ==
=== º Función potencial ===
La velocidad de las partículas del fluido viene dada por el gradiente de la función potencial:
<math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math>
Las curvas de nivel de esta función potencial son las siguientes:
[[Archivo:2curvasdenivel.jpg|800px|thumb|center|FUNCIÓN POTENCIAL]]
El código matlab empleado es:
{{matlab|codigo=
ro=linspace(1,5,50);
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt);
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));
f=(Mro+1./Mro).*sin(Mtt); %función potencial
subplot(1,1,1)
contour(Mx,My,f,50); colorbar %curvas de nivel función potencial
}}

=== º Campo de velocidades (gradiente de función potencial) ===
Si dibujamos nuestro campo de velocidades <math>\vec u=\nabla\varphi</math> en matlab, la representación del mismo con flechas deberá ser ortogonal a las curvas de nivel antes dibujadas. Ya que, de forma geométrica, el gradiente es un vector normal a la curva de nivel de la función de la cual es gradiente en el punto que se esté estudiando.
Hemos pasado nuestra función potencial a cartesianas y trabajado posteriormente con sus derivadas. 
Así, con las igualdades: <math>\rho=\sqrt{x^2+y^2}\qquad\theta=arctg(\frac{y}{x})</math> 
Obtenemos la función potencial: <math>\varphi=y+\frac{y}{x^2+y^2}</math> 

El gradiente se define como:
<math>\vec u=\nabla\varphi=\frac{∂φ}{∂x}\vec i + \frac{∂φ}{∂j}\vec j </math>

Realizamos las derivadas parciales: 
<math>FX=\frac{\partial\varphi(x,y)}{\partial x}=\frac{-2xy}{(x^2+y^2)^2}\qquad FY=\frac{\partial\varphi(x,y)}{\partial y}=\frac{x^2-y^2}{(x^2+y^2)^2}</math> 
Con estas funciones, podremos dibujar en matlab el campo de velocidades con el siguiente código:
{{matlab|codigo=
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));
f=(Mro+1./Mro).*sin(Mtt); %función potencial
syms x y %nuevas variables en cartesianas para definir derivadas (campo de velocidades)
fx=inline('-2*(x.*y)./((x.^2+y.^2).^2)','x','y'); %derivadas que definen nuestro campos de velocidades
fy=inline('1+((x.^2-y.^2)./((x.^2+y.^2).^2))','x','y');
d1=fx(Mx,My);
d2=fy(Mx,My);
hold on
subplot(1,1,1)
contour(Mx,My,f,20); colorbar %curvas de nivel función potencial
quiver(Mx,My,d1,d2) %campo de velocidades fluido, ortogonal a curvas de nivel
axis([-5,5,-5,5])
hold off
}}
En nuestro dibujo saldrán las flechas, que representan nuestro campo de velocidades superpuestas con las curvas de nivel de la función potencial. Si ampliamos la imagen podremos observar que estos dos elementos son ortogonales:

[[Archivo: 2campoycurvas.jpg|800px|thumb|center|CAMPO DE VELOCIDADES]]

==º Ortogonalidad <math>\vec u</math> y <math>\vec n</math> ==
Como el obstáculo es un círculo situado en el plano xy, sabemos con certeza que entonces el vector normal será el correspondiente al eje z, es decir: <math>\vec n= \vec k</math>
Del apartado anterior conocemos que <math>\vec u=\nabla\varphi</math>.

<math>\vec{u}=sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}</math>

Conociendo entonces el vector normal y el vector gradiente, podemos calcular su producto escalar. Para comprobar si dos vectores son ortogonales, su producto escalar debe de ser nulo, es por ello que realizamos el producto escalar del gradiente y del vector normal.
<math>\nabla\varphi\cdot\vec n</math> = <math>(sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}})\cdot\vec ({e}_{z})</math> = 0

Como podemos observar, el producto escalar resulta 0, y queda comprobado así que el vector función potencial es perpendicular a su vector normal.

==º Calculo de <math>\vec u</math> con <math> \frac{1}{\rho}</math> despreciable ==
Al centrarnos en el fluido a una distancia muy lejana, podemos considerar <math> \frac{1}{\rho}</math> despreciable. En este caso sabremos que <math> u</math> cambia. La velocidad de las partículas desde los puntos lejanos del fluido viene definida por el gradiente de la función potencial <math>\varphi=\rho\sin(\theta)</math>

<math>\vec{u}=sin(\theta){\vec{e}_{\rho}}+\frac{1}{\rho}\rho\cos(\theta){\vec{e}_{\theta}}</math> 

En este caso:
<math>|\vec{u}| =\sqrt{cos(\theta)^2+sin(\theta)^2}=1</math>

== º Rotacional y Divergencia ==
:
<math>\nabla\times\vec u=\begin{vmatrix} \vec {e}_{\rho}&\rho\vec {e}_{\theta}&\vec {e}_{z} \\ \frac{\partial}{\partial{\rho}} & \frac{\partial}{\partial{\theta}} & \frac{\partial}{\partial{z}} \\ sin(\theta)[1-\frac{1}{\rho^2}] & \rho\cos(\theta)[\rho +\frac{1}{\rho}] & {0} \end{vmatrix}=[-cos(\theta)(1-\frac{1}{\rho^2})]\vec {e}_{z}+[cos(\theta)(1-\frac{1}{\rho^2})]\vec {e}_{z}=\vec {0}</math>
:
<math>\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}({\rho}sin(\theta)(1 -\frac{1}{\rho^2}))+\frac{\partial}{\partial{\theta}}(cos(\theta)(\rho +\frac{1}{\rho})+\frac{\partial}{\partial{z}}(\rho\cdot{0})]=\frac{1}{\rho}[sin(\theta)(1+\frac{1}{\rho^2})-sin(\theta)(1+\frac{1}{\rho^2})]=0 </math>

En este apartado comprobamos analíticamente que la divergencia y rotacional de <math> u</math> son nulas. Que la divergencia sea nula indica que en este caso el volumen del fluido no varía, es decir, se mantiene, ni se expande ni se contrae. A esta condición se le denomina incompresibilidad. Además, que el rotacional sea nulo demuestra que las partículas del fluido no giran.

== ºLíneas de corriente ==
Las líneas de corriente son aquellas que son tangentes a cada punto del campo <math>\vec u</math>. Para ello calcularemos el vector <math>\vec v=\vec k\times\vec u</math>

<math>\vec v=\frac{1}{\rho}\begin{vmatrix} \vec {e}_{\rho}&\rho\vec {e}_{\theta}&\vec {e}_{z} \\ {0} & {0} & {1} \\ sin(\theta)(1-\frac{1}{\rho^2}) & cos(\theta)(1+\frac{1}{\rho^2})\rho & {0} \end{vmatrix}= sin(\theta)(1-\frac{1}{\rho^2})\vec {e}_{\theta} - cos(\theta)(1+\frac{1}{\rho^2})\vec {e}_{\rho} =\vec {grad}\psi</math>:

<math>\psi=\int - cos(\theta)(1+\frac{1}{\rho^2}) \, d\rho \,\!=-cos(\theta)(\rho-\frac{1}{\rho}) + f(\theta) </math>:

<math>\frac{\partial\psi}{\partial\theta}=sin(\theta)(\rho-\frac{1}{\rho})+f'(\theta)=sin(\theta)(\rho-\frac{1}{\rho})\qquad f'(\theta)=0\qquad f(\theta)=cte</math>:

<math>\psi=-cos(\theta)(\rho-\frac{1}{\rho})</math>:

Del apartado anterior, podemos ver que el rotacional es nulo, y por lo tanto es irrotacional. De esta forma, las líneas de corriente serán tangentes al campo. 
Las curvas de nivel de la corriente son las siguientes: 
[[Archivo:curvascorriente.jpg|800px|miniaturadeimagen|centro|LÍNEAS DE CORRIENTE]] 
Si superponemos la imagen con el campo de velocidades antes calculado, observamos que efectivamente son tangentes: 
[[Archivo: curvasycampo.png|800px|miniaturadeimagen|centro|Representación de las líneas de corriente, tangentes al campo <math>\vec u</math>]]

El código Matlab empleado es el siguiente:

{{matlab|codigo=ro=linspace(1,5,300); %intervalo de ro (1,5)
tt=linspace(0,2*pi,300); %intervalo de teta (0,2pi)
[Mro,Mtt]=meshgrid(ro,tt);
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));

syms x y
fi=inline('-x+x./(x.^2+y.^2)','x','y'); %funciÃ³n que define la corriente de u
f2=fi(Mx,My);
fx=inline('-2*(x.*y)./((x.^2+y.^2).^2)','x','y'); %derivadas de funcion potencial que definen campo de velocidades del fluido
fy=inline('1+((x.^2-y.^2)./((x.^2+y.^2).^2))','x','y');
d1=fx(Mx,My);
d2=fy(Mx,My);
hold on
subplot(1,1,1)
contour(Mx,My,f2,100) %curvas de nivel de la corriente
quiver(Mx,My,d1,d2) %campo de velocidades fluido, tangente a lineas de corriente
hold off
}}

== º Puntos frontera==
Al tener definida la superficie del anillo con los radios 1 y 5, buscamos donde están situados los puntos de frontera. Con la función potencial inicial <math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math> y con
<math> u</math> = <math>(sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}) </math>.

Calculamos el modulo del vector en función de <math>\theta</math> para

<math>\rho=1\qquad </math> <math> u</math> = <math> 2cos(\theta){\vec{e}_{\theta}} </math>

<math>|\vec u|=\sqrt{4cos(\theta)^2} = 2cos(\theta) = 2|cos(\theta)|</math> 

Como podemos observar, el módulo de la velocidad es máximo cuando cos<math>\theta</math> alcanza sus valores máximos (en 0 y en <math>\pi</math>, donde alcanza en valor absoluto el valor de 1). 

Así, los puntos de velocidad máxima son <math>\theta=0</math> y <math>\theta={\pi}</math>. 

Los puntos de velocidad mínima o puntos de remanso son <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math>. 

== º Presión del fluido ==

En dinámica de fluidos, el principio de Bernoulli, también denominado ecuación de Bernoulli, describe el comportamiento de un fluido moviéndose a lo largo de una línea de corriente. Relacionando las siguientes variables: 

<math>p =</math>presión estática a la que está sometido el fluido debida a las moléculas que lo rodean. 
<math>d =</math>densidad del fluido. 

<math>\vec { u } =</math>velocidad de flujo del fluido. :

<math>\frac { 1 }{ 2 } d { \left| \vec { u } \right| }^{ 2 }+ { p } =cte</math> 
Para que esta expresión pueda aplicarse, el fluido debe ser incompresible, no viscoso y fluir en régimen estacionario (la velocidad en un punto determinado no varía con el tiempo). Tomaremos nuestro fluido como tal. 
Suponiendo que nuestro fluido efectivamente cumple la ecuación de Bernouilli, que posee una densidad constante de d=2 y tomando la cte de la ecuación como cte=10, calcularemos la presión de nuestro fluido. 
Sustituyendo los valores mencionados tenemos que: 
<math> p = 10- |\vec{u}|^2 </math> 
Calculamos el módulo al cuadrado de nuestro campo de velocidades: 
<math> |\vec{u}|^2 = (1-\frac{1}{ρ^2})^2 sen^2(θ) + (1+\frac{1}{ρ^2})^2 cos^2(θ) </math> 
Así la función que define la presión del fluido es: 
<math> p = 10- (1-\frac{1}{ρ^2})^2 sen^2(θ) + (1+\frac{1}{ρ^2})^2 cos^2(θ) </math> 

Para la representación de la presión empleamos el siguiente código:
{{matlab| codigo=
ro=linspace(1,5,50);
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt);
x=Mro.*cos(Mtt);
y=Mro.*sin(Mtt);
P=10-(((1-1./Mro.^2).^2).*(sin(Mtt).^2)+((1+1./Mro.^2).^2).*(cos(Mtt).^2)); %función presión
surf(x,y,P); colorbar;
axis([-5,5,-5,5])
view(2)
Pmax=max(max(P)) %presión máx y mín
Pmin=min(min(P))
}}
Matlab nos devuelve la siguiente imagen y los valores máximos y mínimos de la presión (Presión máx: 9.9959 Presión mínima: 6)
[[Archivo:Presionnn.jpg|800px|miniaturadeimagen|centro|Fig.PRESIÓN DEL FLUIDO]] 
Recordemos ahora la imagen del campo de velocidades (ahora sin ampliar) que habíamos obtenido anteriormente con el código matlab de el apartado 2.2: 

[[Archivo:velocidadd.jpg|800px|miniaturadeimagen|centro|Fig.CAMPO DE VELOCIDADES]] 

Si comparamos nuestra imagen de la presión del fluido (Fig.PRESIÓN DEL FLUIDO) con la del campo de velocidades (Fig.CAMPO DE VELOCIDADES) podemos observar que los puntos que cuentan con módulo de velocidad máxima (<math>\theta=0</math> y <math>\theta={\pi}</math>) son aquellos con presiones mínimas. Esto es un resultado coherente que verifica la ecuación de Bernoulli, a mayor velocidad, menores presiones.

== ºEcuación de Navier-Stokes ==

La ecuación de Navier-Stokes describe el movimiento de los fluidos a partir de su campo de velocidades y de presiones:

<math> d(\vec{u}·\nabla)\vec{u} + \nabla p = \mu \Delta \vec{u} </math> 

Los fluidos incompresibles tienen la propiedad de ocupar siempre el mismo volumen, de esta forma la viscosidad será nula: <math> \mu =0 </math> y la densidad será nula también: <math>d = 2</math>. Si sustituimos en la ecuación general, llegamos a la conclusión de que para un fluido incompresible:

<math> 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 </math> 

Si consideramos los campos <math> (\vec{u},p) </math> hallados anteriormente:

<math> (\vec{u}·\nabla)\vec{u} = u_{1} \frac{\partial u_{1}}{\partial \rho} \vec{e}_{\rho} + u_{1} \frac{\partial u_{1}}{\partial \theta} \vec{e}_{\theta}+ u_{2} \frac{\partial u_{2}}{\partial \rho} \vec{e}_{\rho}+ u_{2} \frac{\partial u_{2}}{\partial \theta} \vec{e}_{\theta} = </math> :

<math> (\vec{u}·\nabla)\vec{u}=(sen(\theta)(1-\frac{1}{\rho^2})[\frac{2}{\rho^3}sin(\theta)\vec{e}_{\rho}+(1-\frac{1}{\rho^2})cos(\theta)\vec{e}_{\theta}])+((1+\frac{1}{\rho^2})cos(\theta)[-\frac{2}{\rho^3}cos(\theta)\vec{e}_{\rho}-(1+\frac{1}{\rho^2})sin(\theta)\vec{e}_{\theta}])=</math>

<math> =(\vec{u}·\nabla)\vec{u}=[\frac{2}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)-\frac{2}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[(1-\frac{1}{\rho^2})^2sin(\theta)cos(\theta)-(1+\frac{1}{\rho^2})^2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

Lo aplicamos a nuestro caso con la densidad determinada:

<math> 2(\vec{u}·\nabla)\vec{u}=[\frac{4}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)- \frac{4}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[(1-\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)-(1+\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

De la misma forma:

<math> \nabla p = \frac{\partial p}{\partial \rho} \vec{e}_{\rho} + \frac{\partial p}{\partial \theta} \vec{e}_{\theta} = </math> :

<math> \nabla p = -[(2(1-\frac{1}{\rho^2})\frac{2}{\rho^3}sin^2(\theta))+(2(1+\frac{1}{\rho^2})(-\frac{2}{\rho^3})cos(\theta))]\vec{e}_{\rho} - [((1-\frac{1}{\rho^2})^2)2sin(\theta)cos(\theta)+((1+\frac{1}{\rho^2})^2)2(-sin(\theta))cos(\theta)]\vec{e}_{\theta}=</math>

<math> = [-\frac{4}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)+ \frac{4}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[-(1-\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)+(1+\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

Habiendo sustituido todos los valores en la ecuación de Navier-Stokes, coparando los resultados vemos directamente que se cumple <math> 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 </math> quedando así demostrado que se cumple la ecuación.

== º Trayectoria de la línea de corriente==
Si fuéramos una partícula del fluido con el que estamos trabajando, seguiríamos la trayectoria de una línea de corriente, es decir empezaríamos con su trayectoria original y a medida que se va acercando al obstáculo recorreríamos el contorno de este, y cuando lo hubiésemos pasado, volveríamos a la dirección inicial. En nuestro caso, como el obstáculo tiene un contorno circular el fluido empieza recto, luego se curva para evitar el obstáculo y finalmente vuelve a ir en la dirección inicial.
En relación a la variación de la velocidad y presión podemos comprobar con las gráficas, que en los puntos donde la velocidad es mínima las presiones son máximas. En los puntos que ya hemos calculado anteriormente donde la velocidad es 0, <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math>, sabemos que la presión es máxima.
Con esto podemos deducir que la presión y velocidad de estas partículas son inversamente proporcionales.

== º Paradoja D'Alembert==
:
<math>\vec{u}=sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}</math>

Para <math>\rho=1</math>:

<math>\vec u=2cos(\theta)\vec {e}_{\theta}</math>:

Según el Teorema de Kutta-Joukowski la fuerza que ejerce el fluido sobre el obstáculo, en este caso el círculo de radio 1, es proporcional a la circulación. Por lo tanto:

<math>\int_{0}^{2\pi} \vec u r'(\theta)\,d\theta =\int_{0}^{2\pi} 2cos(\theta)\vec{e}_{\theta}\vec{e}_{\theta}\,d\theta =\int_{0}^{2\pi} 2cos(\theta)\,d\theta=2[sin(0)-sin(2\pi)]=0 </math>:

El fluido no ejerce ninguna fuerza sobre el obstáculo, dando lugar a la paradoja D'Alembert, que concluye que la fuerza resultante sobre el cuerpo es cero, contradiciéndose con la observación. Esto demuestra también que la circulación del campo vectorial a lo largo de la circunferencia se anula.
Para poder ilustrar esta fuerza se necesitan tomar soluciones de modelos de fluidos con una viscosidad positiva.

== º Curvas de nivel de la presión==
Observando las curvas de nivel de la presión podemos entender mejor las conclusiones a las que hemos ido llegando a lo largo del trabajo. En esta gráfica podemos observar los puntos de presión máxima <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math> y los puntos de presión mínima <math>\theta=0</math> y <math>\theta=\pi</math>.
Además, cabe destacar que, en las zonas en las que las curvas de nivel están mas juntas, hay mayor variación de la presión.

[[Archivo:Doce_mnj.jpg|800px|miniaturadeimagen|centro|Líneas de presión]]

Código Matlab:
{{matlab| codigo=

contour(x,y,P,100) % Curvas de nivel de la presión
axis([-5,5,-5,5]) % Región a dibujar

}}
[[Categoría:Teoría de Campos]]
[[Categoría:TC21/22]]

Visualización de campos escalares y vectoriales en fluidos. (Grupo B6)

2021-12-13T07:17:20Z

Julia Cabeza:

{{ TrabajoED | Visualización de campos escalares y vectoriales en fluidos. (Grupo B6) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC21/22|2021-22]] | María Magaldi Jurado, Nuria Martínez Ballester, Julia Cabeza Duque }}

En este trabajo vamos a visualizar campos escalares y vectoriales en fluidos. Más específicamente, el fluido a tratar es uno incompresible que se encuentra alrededor de un obstáculo de forma circular, formando así un anillo comprendido entre las circunferencias de radio 1 y 5, con el centro en el origen de coordenadas. Para facilitar el trabajo usaremos coordenadas polares (cilíndricas <math>(\rho,\theta)</math>). Para dibujar que el fluido va por el exterior de la superficie definida consideraremos que <math>(\rho,\theta) \in [0,5]*[0,2\pi]</math>.
Además, tendremos que la velocidad de las partículas del fluido viene definida por el gradiente de la función potencial <math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math>

== º Mallado del fluido ==
El siguiente mallado representa los puntos interiores de la región ocupada por un fluido. Se encuentra en el exterior del círculo unidad. Como ya se ha introducido, nuestro mallado del anillo está comprendido entre los radios 1 y 5 y centro el origen.

[[Archivo:1,mallad0.jpg|800px|thumb|center|Mallado del fluido]]

El CÓDIGO matlab empleado es:
{{matlab|codigo=
ro=linspace(1,5,50); %creación de los intervalos
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt); %creación matrices
Mx=Mro.*cos(Mtt); %parametrización de x
My=Mro.*sin(Mtt); %parametrización de y
Mz=zeros(size(Mx)); %parametrización de z (nula al estar en un plano)
subplot(1,1,1);
mesh(Mx,My,Mz) %mallado
axis([-5,5,-5,5]) %superficie definida en [-5,5]*[-5,5]
view(2)
}}

== º Función potencial <math>\varphi</math> y velocidad de las partículas del fluido ==
=== º Función potencial ===
La velocidad de las partículas del fluido viene dada por el gradiente de la función potencial:
<math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math>
Las curvas de nivel de esta función potencial son las siguientes:
[[Archivo:2curvasdenivel.jpg|800px|thumb|center|FUNCIÓN POTENCIAL]]
El código matlab empleado es:
{{matlab|codigo=
ro=linspace(1,5,50);
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt);
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));
f=(Mro+1./Mro).*sin(Mtt); %función potencial
subplot(1,1,1)
contour(Mx,My,f,50); colorbar %curvas de nivel función potencial
}}

=== º Campo de velocidades (gradiente de función potencial) ===
Si dibujamos nuestro campo de velocidades <math>\vec u=\nabla\varphi</math> en matlab, la representación del mismo con flechas deberá ser ortogonal a las curvas de nivel antes dibujadas. Ya que, de forma geométrica, el gradiente es un vector normal a la curva de nivel de la función de la cual es gradiente en el punto que se esté estudiando.
Hemos pasado nuestra función potencial a cartesianas y trabajado posteriormente con sus derivadas. 
Así, con las igualdades: <math>\rho=\sqrt{x^2+y^2}\qquad\theta=arctg(\frac{y}{x})</math> 
Obtenemos la función potencial: <math>\varphi=y+\frac{y}{x^2+y^2}</math> 

El gradiente se define como:
<math>\vec u=\nabla\varphi=\frac{∂φ}{∂x}\vec i + \frac{∂φ}{∂j}\vec j </math>

Realizamos las derivadas parciales: 
<math>FX=\frac{\partial\varphi(x,y)}{\partial x}=\frac{-2xy}{(x^2+y^2)^2}\qquad FY=\frac{\partial\varphi(x,y)}{\partial y}=\frac{x^2-y^2}{(x^2+y^2)^2}</math> 
Con estas funciones, podremos dibujar en matlab el campo de velocidades con el siguiente código:
{{matlab|codigo=
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));
f=(Mro+1./Mro).*sin(Mtt); %función potencial
syms x y %nuevas variables en cartesianas para definir derivadas (campo de velocidades)
fx=inline('-2*(x.*y)./((x.^2+y.^2).^2)','x','y'); %derivadas que definen nuestro campos de velocidades
fy=inline('1+((x.^2-y.^2)./((x.^2+y.^2).^2))','x','y');
d1=fx(Mx,My);
d2=fy(Mx,My);
hold on
subplot(1,1,1)
contour(Mx,My,f,20); colorbar %curvas de nivel función potencial
quiver(Mx,My,d1,d2) %campo de velocidades fluido, ortogonal a curvas de nivel
axis([-5,5,-5,5])
hold off
}}
En nuestro dibujo saldrán las flechas, que representan nuestro campo de velocidades superpuestas con las curvas de nivel de la función potencial. Si ampliamos la imagen podremos observar que estos dos elementos son ortogonales:

[[Archivo: 2campoycurvas.jpg|800px|thumb|center|CAMPO DE VELOCIDADES]]

==º Ortogonalidad <math>\vec u</math> y <math>\vec n</math> ==
Como el obstáculo es un círculo situado en el plano xy, sabemos con certeza que entonces el vector normal será el correspondiente al eje z, es decir: <math>\vec n= \vec k</math>
Del apartado anterior conocemos que <math>\vec u=\nabla\varphi</math>.

<math>\vec{u}=sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}</math>

Conociendo entonces el vector normal y el vector gradiente, podemos calcular su producto escalar. Para comprobar si dos vectores son ortogonales, su producto escalar debe de ser nulo, es por ello que realizamos el producto escalar del gradiente y del vector normal.
<math>\nabla\varphi\cdot\vec n</math> = <math>(sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}})\cdot\vec ({e}_{z})</math> = 0

Como podemos observar, el producto escalar resulta 0, y queda comprobado así que el vector función potencial es perpendicular a su vector normal.

==º Calculo de <math>\vec u</math> con <math> \frac{1}{\rho}</math> despreciable ==
Al centrarnos en el fluido a una distancia muy lejana, podemos considerar <math> \frac{1}{\rho}</math> despreciable. En este caso sabremos que <math> u</math> cambia. La velocidad de las partículas desde los puntos lejanos del fluido viene definida por el gradiente de la función potencial <math>\varphi=\rho\sin(\theta)</math>

<math>\vec{u}=sin(\theta){\vec{e}_{\rho}}+\frac{1}{\rho}\rho\cos(\theta){\vec{e}_{\theta}}</math> 

En este caso:
<math>|\vec{u}| =\sqrt{cos(\theta)^2+sin(\theta)^2}=1</math>

== º Rotacional y Divergencia ==
:
<math>\nabla\times\vec u=\begin{vmatrix} \vec {e}_{\rho}&\rho\vec {e}_{\theta}&\vec {e}_{z} \\ \frac{\partial}{\partial{\rho}} & \frac{\partial}{\partial{\theta}} & \frac{\partial}{\partial{z}} \\ sin(\theta)[1-\frac{1}{\rho^2}] & \rho\cos(\theta)[\rho +\frac{1}{\rho}] & {0} \end{vmatrix}=[-cos(\theta)(1-\frac{1}{\rho^2})]\vec {e}_{z}+[cos(\theta)(1-\frac{1}{\rho^2})]\vec {e}_{z}=\vec {0}</math>
:
<math>\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}({\rho}sin(\theta)(1 -\frac{1}{\rho^2}))+\frac{\partial}{\partial{\theta}}(cos(\theta)(\rho +\frac{1}{\rho})+\frac{\partial}{\partial{z}}(\rho\cdot{0})]=\frac{1}{\rho}[sin(\theta)(1+\frac{1}{\rho^2})-sin(\theta)(1+\frac{1}{\rho^2})]=0 </math>

En este apartado comprobamos analíticamente que la divergencia y rotacional de <math> u</math> son nulas. Que la divergencia sea nula indica que en este caso el volumen del fluido no varía, es decir, se mantiene, ni se expande ni se contrae. A esta condición se le denomina incompresibilidad. Además, que el rotacional sea nulo demuestra que las partículas del fluido no giran.

== ºLíneas de corriente ==
Las líneas de corriente son aquellas que son tangentes a cada punto del campo <math>\vec u</math>. Para ello calcularemos el vector <math>\vec v=\vec k\times\vec u</math>

<math>\vec v=\frac{1}{\rho}\begin{vmatrix} \vec {e}_{\rho}&\rho\vec {e}_{\theta}&\vec {e}_{z} \\ {0} & {0} & {1} \\ sin(\theta)(1-\frac{1}{\rho^2}) & cos(\theta)(1+\frac{1}{\rho^2})\rho & {0} \end{vmatrix}= sin(\theta)(1-\frac{1}{\rho^2})\vec {e}_{\theta} - cos(\theta)(1+\frac{1}{\rho^2})\vec {e}_{\rho} =\vec {grad}\psi</math>:

<math>\psi=\int - cos(\theta)(1+\frac{1}{\rho^2}) \, d\rho \,\!=-cos(\theta)(\rho-\frac{1}{\rho}) + f(\theta) </math>:

<math>\frac{\partial\psi}{\partial\theta}=sin(\theta)(\rho-\frac{1}{\rho})+f'(\theta)=sin(\theta)(\rho-\frac{1}{\rho})\qquad f'(\theta)=0\qquad f(\theta)=cte</math>:

<math>\psi=-cos(\theta)(\rho-\frac{1}{\rho})</math>:

Del apartado anterior, podemos ver que el rotacional es nulo, y por lo tanto es irrotacional. De esta forma, las líneas de corriente serán tangentes al campo. 
Las curvas de nivel de la corriente son las siguientes: 
[[Archivo:curvascorriente.jpg|800px|miniaturadeimagen|centro|LÍNEAS DE CORRIENTE]] 
Si superponemos la imagen con el campo de velocidades antes calculado, observamos que efectivamente son tangentes: 
[[Archivo: curvasycampo.png|800px|miniaturadeimagen|centro|Representación de las líneas de corriente, tangentes al campo <math>\vec u</math>]]

El código Matlab empleado es el siguiente:

{{matlab|codigo=ro=linspace(1,5,300); %intervalo de ro (1,5)
tt=linspace(0,2*pi,300); %intervalo de teta (0,2pi)
[Mro,Mtt]=meshgrid(ro,tt);
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));

syms x y
fi=inline('-x+x./(x.^2+y.^2)','x','y'); %funciÃ³n que define la corriente de u
f2=fi(Mx,My);
fx=inline('-2*(x.*y)./((x.^2+y.^2).^2)','x','y'); %derivadas de funcion potencial que definen campo de velocidades del fluido
fy=inline('1+((x.^2-y.^2)./((x.^2+y.^2).^2))','x','y');
d1=fx(Mx,My);
d2=fy(Mx,My);
hold on
subplot(1,1,1)
contour(Mx,My,f2,100) %curvas de nivel de la corriente
quiver(Mx,My,d1,d2) %campo de velocidades fluido, tangente a lineas de corriente
hold off
}}

== º Puntos frontera==
Al tener definida la superficie del anillo con los radios 1 y 5, buscamos donde están situados los puntos de frontera. Con la función potencial inicial <math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math> y con
<math> u</math> = <math>(sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}) </math>.

Calculamos el modulo del vector en función de <math>\theta</math> para

<math>\rho=1\qquad </math> <math> u</math> = <math> 2cos(\theta){\vec{e}_{\theta}} </math>

<math>|\vec u|=\sqrt{4cos(\theta)^2} = 2cos(\theta) = 2|cos(\theta)|</math> 

Como podemos observar, el módulo de la velocidad es máximo cuando cos<math>\theta</math> alcanza sus valores máximos (en 0 y en <math>\pi</math>, donde alcanza en valor absoluto el valor de 1). 

Así, los puntos de velocidad máxima son <math>\theta=0</math> y <math>\theta={pi}</math>. 

Los puntos de velocidad mínima o puntos de remanso son <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math>. 

== º Presión del fluido ==

En dinámica de fluidos, el principio de Bernoulli, también denominado ecuación de Bernoulli, describe el comportamiento de un fluido moviéndose a lo largo de una línea de corriente. Relacionando las siguientes variables: 

<math>p =</math>presión estática a la que está sometido el fluido debida a las moléculas que lo rodean. 
<math>d =</math>densidad del fluido. 

<math>\vec { u } =</math>velocidad de flujo del fluido. :

<math>\frac { 1 }{ 2 } d { \left| \vec { u } \right| }^{ 2 }+ { p } =cte</math> 
Para que esta expresión pueda aplicarse, el fluido debe ser incompresible, no viscoso y fluir en régimen estacionario (la velocidad en un punto determinado no varía con el tiempo). Tomaremos nuestro fluido como tal. 
Suponiendo que nuestro fluido efectivamente cumple la ecuación de Bernouilli, que posee una densidad constante de d=2 y tomando la cte de la ecuación como cte=10, calcularemos la presión de nuestro fluido. 
Sustituyendo los valores mencionados tenemos que: 
<math> p = 10- |\vec{u}|^2 </math> 
Calculamos el módulo al cuadrado de nuestro campo de velocidades: 
<math> |\vec{u}|^2 = (1-\frac{1}{ρ^2})^2 sen^2(θ) + (1+\frac{1}{ρ^2})^2 cos^2(θ) </math> 
Así la función que define la presión del fluido es: 
<math> p = 10- (1-\frac{1}{ρ^2})^2 sen^2(θ) + (1+\frac{1}{ρ^2})^2 cos^2(θ) </math> 

Para la representación de la presión empleamos el siguiente código:
{{matlab| codigo=
ro=linspace(1,5,50);
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt);
x=Mro.*cos(Mtt);
y=Mro.*sin(Mtt);
P=10-(((1-1./Mro.^2).^2).*(sin(Mtt).^2)+((1+1./Mro.^2).^2).*(cos(Mtt).^2)); %función presión
surf(x,y,P); colorbar;
axis([-5,5,-5,5])
view(2)
Pmax=max(max(P)) %presión máx y mín
Pmin=min(min(P))
}}
Matlab nos devuelve la siguiente imagen y los valores máximos y mínimos de la presión (Presión máx: 9.9959 Presión mínima: 6)
[[Archivo:Presionnn.jpg|800px|miniaturadeimagen|centro|Fig.PRESIÓN DEL FLUIDO]] 
Recordemos ahora la imagen del campo de velocidades (ahora sin ampliar) que habíamos obtenido anteriormente con el código matlab de el apartado 2.2: 

[[Archivo:velocidadd.jpg|800px|miniaturadeimagen|centro|Fig.CAMPO DE VELOCIDADES]] 

Si comparamos nuestra imagen de la presión del fluido (Fig.PRESIÓN DEL FLUIDO) con la del campo de velocidades (Fig.CAMPO DE VELOCIDADES) podemos observar que los puntos que cuentan con módulo de velocidad máxima (<math>\theta=0</math> y <math>\theta={\pi}</math>) son aquellos con presiones mínimas. Esto es un resultado coherente que verifica la ecuación de Bernoulli, a mayor velocidad, menores presiones.

== ºEcuación de Navier-Stokes ==

La ecuación de Navier-Stokes describe el movimiento de los fluidos a partir de su campo de velocidades y de presiones:

<math> d(\vec{u}·\nabla)\vec{u} + \nabla p = \mu \Delta \vec{u} </math> 

Los fluidos incompresibles tienen la propiedad de ocupar siempre el mismo volumen, de esta forma la viscosidad será nula: <math> \mu =0 </math> y la densidad será nula también: <math>d = 2</math>. Si sustituimos en la ecuación general, llegamos a la conclusión de que para un fluido incompresible:

<math> 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 </math> 

Si consideramos los campos <math> (\vec{u},p) </math> hallados anteriormente:

<math> (\vec{u}·\nabla)\vec{u} = u_{1} \frac{\partial u_{1}}{\partial \rho} \vec{e}_{\rho} + u_{1} \frac{\partial u_{1}}{\partial \theta} \vec{e}_{\theta}+ u_{2} \frac{\partial u_{2}}{\partial \rho} \vec{e}_{\rho}+ u_{2} \frac{\partial u_{2}}{\partial \theta} \vec{e}_{\theta} = </math> :

<math> (\vec{u}·\nabla)\vec{u}=(sen(\theta)(1-\frac{1}{\rho^2})[\frac{2}{\rho^3}sin(\theta)\vec{e}_{\rho}+(1-\frac{1}{\rho^2})cos(\theta)\vec{e}_{\theta}])+((1+\frac{1}{\rho^2})cos(\theta)[-\frac{2}{\rho^3}cos(\theta)\vec{e}_{\rho}-(1+\frac{1}{\rho^2})sin(\theta)\vec{e}_{\theta}])=</math>

<math> =(\vec{u}·\nabla)\vec{u}=[\frac{2}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)-\frac{2}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[(1-\frac{1}{\rho^2})^2sin(\theta)cos(\theta)-(1+\frac{1}{\rho^2})^2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

Lo aplicamos a nuestro caso con la densidad determinada:

<math> 2(\vec{u}·\nabla)\vec{u}=[\frac{4}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)- \frac{4}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[(1-\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)-(1+\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

De la misma forma:

<math> \nabla p = \frac{\partial p}{\partial \rho} \vec{e}_{\rho} + \frac{\partial p}{\partial \theta} \vec{e}_{\theta} = </math> :

<math> \nabla p = -[(2(1-\frac{1}{\rho^2})\frac{2}{\rho^3}sin^2(\theta))+(2(1+\frac{1}{\rho^2})(-\frac{2}{\rho^3})cos(\theta))]\vec{e}_{\rho} - [((1-\frac{1}{\rho^2})^2)2sin(\theta)cos(\theta)+((1+\frac{1}{\rho^2})^2)2(-sin(\theta))cos(\theta)]\vec{e}_{\theta}=</math>

<math> = [-\frac{4}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)+ \frac{4}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[-(1-\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)+(1+\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

Habiendo sustituido todos los valores en la ecuación de Navier-Stokes, coparando los resultados vemos directamente que se cumple <math> 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 </math> quedando así demostrado que se cumple la ecuación.

== º Trayectoria de la línea de corriente==
Si fuéramos una partícula del fluido con el que estamos trabajando, seguiríamos la trayectoria de una línea de corriente, es decir empezaríamos con su trayectoria original y a medida que se va acercando al obstáculo recorreríamos el contorno de este, y cuando lo hubiésemos pasado, volveríamos a la dirección inicial. En nuestro caso, como el obstáculo tiene un contorno circular el fluido empieza recto, luego se curva para evitar el obstáculo y finalmente vuelve a ir en la dirección inicial.
En relación a la variación de la velocidad y presión podemos comprobar con las gráficas, que en los puntos donde la velocidad es mínima las presiones son máximas. En los puntos que ya hemos calculado anteriormente donde la velocidad es 0, <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math>, sabemos que la presión es máxima.
Con esto podemos deducir que la presión y velocidad de estas partículas son inversamente proporcionales.

== º Paradoja D'Alembert==
:
<math>\vec{u}=sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}</math>

Para <math>\rho=1</math>:

<math>\vec u=2cos(\theta)\vec {e}_{\theta}</math>:

Según el Teorema de Kutta-Joukowski la fuerza que ejerce el fluido sobre el obstáculo, en este caso el círculo de radio 1, es proporcional a la circulación. Por lo tanto:

<math>\int_{0}^{2\pi} \vec u r'(\theta)\,d\theta =\int_{0}^{2\pi} 2cos(\theta)\vec{e}_{\theta}\vec{e}_{\theta}\,d\theta =\int_{0}^{2\pi} 2cos(\theta)\,d\theta=2[sin(0)-sin(2\pi)]=0 </math>:

El fluido no ejerce ninguna fuerza sobre el obstáculo, dando lugar a la paradoja D'Alembert, que concluye que la fuerza resultante sobre el cuerpo es cero, contradiciéndose con la observación. Esto demuestra también que la circulación del campo vectorial a lo largo de la circunferencia se anula.
Para poder ilustrar esta fuerza se necesitan tomar soluciones de modelos de fluidos con una viscosidad positiva.

== º Curvas de nivel de la presión==
Observando las curvas de nivel de la presión podemos entender mejor las conclusiones a las que hemos ido llegando a lo largo del trabajo. En esta gráfica podemos observar los puntos de presión máxima <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math> y los puntos de presión mínima <math>\theta=0</math> y <math>\theta=\pi</math>.
Además, cabe destacar que, en las zonas en las que las curvas de nivel están mas juntas, hay mayor variación de la presión.

[[Archivo:Doce_mnj.jpg|800px|miniaturadeimagen|centro|Líneas de presión]]

Código Matlab:
{{matlab| codigo=

contour(x,y,P,100) % Curvas de nivel de la presión
axis([-5,5,-5,5]) % Región a dibujar

}}
[[Categoría:Teoría de Campos]]
[[Categoría:TC21/22]]

Visualización de campos escalares y vectoriales en fluidos. (Grupo B6)

2021-12-11T15:50:20Z

Julia Cabeza: /* º Trayectoria de la línea de corriente */

{{ TrabajoED | Visualización de campos escalares y vectoriales en fluidos. (Grupo B6) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC21/22|2021-22]] | María Magaldi Jurado, Nuria Martínez Ballester, Julia Cabeza Duque }}

En este trabajo vamos a visualizar campos escalares y vectoriales en fluidos. Más específicamente, el fluido a tratar es uno incompresible que se encuentra alrededor de un obstáculo de forma circular, formando así un anillo comprendido entre las circunferencias de radio 1 y 5, con el centro en el origen de coordenadas. Para facilitar el trabajo usaremos coordenadas polares (cilíndricas <math>(\rho,\theta)</math>). Para dibujar que el fluido va por el exterior de la superficie definida consideraremos que <math>(\rho,\theta) \in [-5,5]*[-5,5]</math>.
Además, tendremos que la velocidad de las partículas del fluido viene definida por el gradiente de la función potencial <math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math>

== º Mallado del fluido ==
El siguiente mallado representa los puntos interiores de la región ocupada por un fluido. Se encuentra en el exterior del círculo unidad. Como ya se ha introducido, nuestro mallado del anillo está comprendido entre los radios 1 y 5 y centro el origen.

[[Archivo:1,mallad0.jpg|800px|thumb|center|Mallado del fluido]]

El CÓDIGO matlab empleado es:
{{matlab|codigo=
ro=linspace(1,5,50); %creación de los intervalos
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt); %creación matrices
Mx=Mro.*cos(Mtt); %parametrización de x
My=Mro.*sin(Mtt); %parametrización de y
Mz=zeros(size(Mx)); %parametrización de z (nula al estar en un plano)
subplot(1,1,1);
mesh(Mx,My,Mz) %mallado
axis([-5,5,-5,5]) %superficie definida en [-5,5]*[-5,5]
view(2)
}}

== º Función potencial <math>\varphi</math> y velocidad de las partículas del fluido ==
=== º Función potencial ===
La velocidad de las partículas del fluido viene dada por el gradiente de la función potencial:
<math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math>
Las curvas de nivel de esta función potencial son las siguientes:
[[Archivo:2curvasdenivel.jpg|800px|thumb|center|FUNCIÓN POTENCIAL]]
El código matlab empleado es:
{{matlab|codigo=
ro=linspace(1,5,50);
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt);
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));
f=(Mro+1./Mro).*sin(Mtt); %función potencial
subplot(1,1,1)
contour(Mx,My,f,50); colorbar %curvas de nivel función potencial
}}

=== º Campo de velocidades (gradiente de función potencial) ===
Si dibujamos nuestro campo de velocidades <math>\vec u=\nabla\varphi</math> en matlab, la representación del mismo con flechas deberá ser ortogonal a las curvas de nivel antes dibujadas. Ya que, de forma geométrica, el gradiente es un vector normal a la curva de nivel de la función de la cual es gradiente en el punto que se esté estudiando.
Hemos pasado nuestra función potencial a cartesianas y trabajado posteriormente con sus derivadas. 
Así, con las igualdades: <math>\rho=\sqrt{x^2+y^2}\qquad\theta=arctg(\frac{y}{x})</math> 
Obtenemos la función potencial: <math>\varphi=y+\frac{y}{x^2+y^2}</math> 

El gradiente se define como:
<math>\vec u=\nabla\varphi=\frac{∂φ}{∂x}\vec i + \frac{∂φ}{∂j}\vec j </math>

Realizamos las derivadas parciales: 
<math>FX=\frac{\partial\varphi(x,y)}{\partial x}=\frac{-2xy}{(x^2+y^2)^2}\qquad FY=\frac{\partial\varphi(x,y)}{\partial y}=\frac{x^2-y^2}{(x^2+y^2)^2}</math> 
Con estas funciones, podremos dibujar en matlab el campo de velocidades con el siguiente código:
{{matlab|codigo=
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));
f=(Mro+1./Mro).*sin(Mtt); %función potencial
syms x y %nuevas variables en cartesianas para definir derivadas (campo de velocidades)
fx=inline('-2*(x.*y)./((x.^2+y.^2).^2)','x','y'); %derivadas que definen nuestro campos de velocidades
fy=inline('1+((x.^2-y.^2)./((x.^2+y.^2).^2))','x','y');
d1=fx(Mx,My);
d2=fy(Mx,My);
hold on
subplot(1,1,1)
contour(Mx,My,f,20); colorbar %curvas de nivel función potencial
quiver(Mx,My,d1,d2) %campo de velocidades fluido, ortogonal a curvas de nivel
axis([-5,5,-5,5])
hold off
}}
En nuestro dibujo saldrán las flechas, que representan nuestro campo de velocidades superpuestas con las curvas de nivel de la función potencial. Si ampliamos la imagen podremos observar que estos dos elementos son ortogonales:

[[Archivo: 2campoycurvas.jpg|800px|thumb|center|CAMPO DE VELOCIDADES]]

==º Ortogonalidad <math>\vec u</math> y <math>\vec n</math> ==
Como el obstáculo es un círculo situado en el plano xy, sabemos con certeza que entonces el vector normal será el correspondiente al eje z, es decir: <math>\vec n= \vec k</math>
Del apartado anterior conocemos que <math>\vec u=\nabla\varphi</math>.

<math>\vec{u}=sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}</math>

Conociendo entonces el vector normal y el vector gradiente, podemos calcular su producto escalar. Para comprobar si dos vectores son ortogonales, su producto escalar debe de ser nulo, es por ello que realizamos el producto escalar del gradiente y del vector normal.
<math>\nabla\varphi\cdot\vec n</math> = <math>(sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}})\cdot\vec ({e}_{z})</math> = 0

Como podemos observar, el producto escalar resulta 0, y queda comprobado así que el vector función potencial es perpendicular a su vector normal.

==º Calculo de <math>\vec u</math> con <math> \frac{1}{\rho}</math> despreciable ==
Al centrarnos en el fluido a una distancia muy lejana, podemos considerar <math> \frac{1}{\rho}</math> despreciable. En este caso sabremos que <math> u</math> cambia. La velocidad de las partículas desde los puntos lejanos del fluido viene definida por el gradiente de la función potencial <math>\varphi=\rho\sin(\theta)</math>

<math>\vec{u}=sin(\theta){\vec{e}_{\rho}}+\frac{1}{\rho}\rho\cos(\theta){\vec{e}_{\theta}}</math> 

En este caso:
<math>|\vec{u}| =\sqrt{cos(\theta)^2+sin(\theta)^2}=1</math>

== º Rotacional y Divergencia ==
:
<math>\nabla\times\vec u=\begin{vmatrix} \vec {e}_{\rho}&\rho\vec {e}_{\theta}&\vec {e}_{z} \\ \frac{\partial}{\partial{\rho}} & \frac{\partial}{\partial{\theta}} & \frac{\partial}{\partial{z}} \\ sin(\theta)[1-\frac{1}{\rho^2}] & \rho\cos(\theta)[\rho +\frac{1}{\rho}] & {0} \end{vmatrix}=[-cos(\theta)(1-\frac{1}{\rho^2})]\vec {e}_{z}+[cos(\theta)(1-\frac{1}{\rho^2})]\vec {e}_{z}=\vec {0}</math>
:
<math>\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}({\rho}sin(\theta)(1 -\frac{1}{\rho^2}))+\frac{\partial}{\partial{\theta}}(cos(\theta)(\rho +\frac{1}{\rho})+\frac{\partial}{\partial{z}}(\rho\cdot{0})]=\frac{1}{\rho}[sin(\theta)(1+\frac{1}{\rho^2})-sin(\theta)(1+\frac{1}{\rho^2})]=0 </math>

En este apartado comprobamos analíticamente que la divergencia y rotacional de <math> u</math> son nulas. Que la divergencia sea nula indica que en este caso el volumen del fluido no varía, es decir, se mantiene, ni se expande ni se contrae. A esta condición se le denomina incompresibilidad. Además, que el rotacional sea nulo demuestra que las partículas del fluido no giran.

== ºLíneas de corriente ==
Las líneas de corriente son aquellas que son tangentes a cada punto del campo <math>\vec u</math>. Para ello calcularemos el vector <math>\vec v=\vec k\times\vec u</math>

<math>\vec v=\frac{1}{\rho}\begin{vmatrix} \vec {e}_{\rho}&\rho\vec {e}_{\theta}&\vec {e}_{z} \\ {0} & {0} & {1} \\ sin(\theta)(1-\frac{1}{\rho^2}) & cos(\theta)(1+\frac{1}{\rho^2})\rho & {0} \end{vmatrix}= sin(\theta)(1-\frac{1}{\rho^2})\vec {e}_{\theta} - cos(\theta)(1+\frac{1}{\rho^2})\vec {e}_{\rho} =\vec {grad}\psi</math>:

<math>\psi=\int - cos(\theta)(1+\frac{1}{\rho^2}) \, d\rho \,\!=-cos(\theta)(\rho-\frac{1}{\rho}) + f(\theta) </math>:

<math>\frac{\partial\psi}{\partial\theta}=sin(\theta)(\rho-\frac{1}{\rho})+f'(\theta)=sin(\theta)(\rho-\frac{1}{\rho})\qquad f'(\theta)=0\qquad f(\theta)=cte</math>:

<math>\psi=-cos(\theta)(\rho-\frac{1}{\rho})</math>:

Del apartado anterior, podemos ver que el rotacional es nulo, y por lo tanto es irrotacional. De esta forma, las líneas de corriente serán tangentes al campo. 
Las curvas de nivel de la corriente son las siguientes: 
[[Archivo:curvascorriente.jpg|800px|miniaturadeimagen|centro|LÍNEAS DE CORRIENTE]] 
Si superponemos la imagen con el campo de velocidades antes calculado, observamos que efectivamente son tangentes: 
[[Archivo: curvasycampo.png|800px|miniaturadeimagen|centro|Representación de las líneas de corriente, tangentes al campo <math>\vec u</math>]]

El código Matlab empleado es el siguiente:

{{matlab|codigo=ro=linspace(1,5,300); %intervalo de ro (1,5)
tt=linspace(0,2*pi,300); %intervalo de teta (0,2pi)
[Mro,Mtt]=meshgrid(ro,tt);
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));

syms x y
fi=inline('-x+x./(x.^2+y.^2)','x','y'); %funciÃ³n que define la corriente de u
f2=fi(Mx,My);
fx=inline('-2*(x.*y)./((x.^2+y.^2).^2)','x','y'); %derivadas de funcion potencial que definen campo de velocidades del fluido
fy=inline('1+((x.^2-y.^2)./((x.^2+y.^2).^2))','x','y');
d1=fx(Mx,My);
d2=fy(Mx,My);
hold on
subplot(1,1,1)
contour(Mx,My,f2,100) %curvas de nivel de la corriente
quiver(Mx,My,d1,d2) %campo de velocidades fluido, tangente a lineas de corriente
hold off
}}

== º Puntos frontera==
Al tener definida la superficie del anillo con los radios 1 y 5, buscamos donde están situados los puntos de frontera. Con la función potencial inicial <math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math> y con
<math> u</math> = <math>(sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}) </math>.

Calculamos el modulo del vector en función de <math>\theta</math> para

<math>\rho=1\qquad </math> <math> u</math> = <math> 2cos(\theta){\vec{e}_{\theta}} </math>

<math>|\vec u|=\sqrt{4cos(\theta)^2} = 2cos(\theta) = 2|cos(\theta)|</math> 

Como podemos observar, el módulo de la velocidad es máximo cuando cos<math>\theta</math> alcanza sus valores máximos (en 0 y en <math>\pi</math>, donde alcanza en valor absoluto el valor de 1). 

Así, los puntos de velocidad máxima son <math>\theta=0</math> y <math>\theta={pi}</math>. 

Los puntos de velocidad mínima o puntos de remanso son <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math>. 

== º Presión del fluido ==

En dinámica de fluidos, el principio de Bernoulli, también denominado ecuación de Bernoulli, describe el comportamiento de un fluido moviéndose a lo largo de una línea de corriente. Relacionando las siguientes variables: 

<math>p =</math>presión estática a la que está sometido el fluido debida a las moléculas que lo rodean. 
<math>d =</math>densidad del fluido. 

<math>\vec { u } =</math>velocidad de flujo del fluido. :

<math>\frac { 1 }{ 2 } d { \left| \vec { u } \right| }^{ 2 }+ { p } =cte</math> 
Para que esta expresión pueda aplicarse, el fluido debe ser incompresible, no viscoso y fluir en régimen estacionario (la velocidad en un punto determinado no varía con el tiempo). Tomaremos nuestro fluido como tal. 
Suponiendo que nuestro fluido efectivamente cumple la ecuación de Bernouilli, que posee una densidad constante de d=2 y tomando la cte de la ecuación como cte=10, calcularemos la presión de nuestro fluido. 
Sustituyendo los valores mencionados tenemos que: 
<math> p = 10- |\vec{u}|^2 </math> 
Calculamos el módulo al cuadrado de nuestro campo de velocidades: 
<math> |\vec{u}|^2 = (1-\frac{1}{ρ^2})^2 sen^2(θ) + (1+\frac{1}{ρ^2})^2 cos^2(θ) </math> 
Así la función que define la presión del fluido es: 
<math> p = 10- (1-\frac{1}{ρ^2})^2 sen^2(θ) + (1+\frac{1}{ρ^2})^2 cos^2(θ) </math> 

Para la representación de la presión empleamos el siguiente código:
{{matlab| codigo=
ro=linspace(1,5,50);
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt);
x=Mro.*cos(Mtt);
y=Mro.*sin(Mtt);
P=10-(((1-1./Mro.^2).^2).*(sin(Mtt).^2)+((1+1./Mro.^2).^2).*(cos(Mtt).^2)); %función presión
surf(x,y,P); colorbar;
axis([-5,5,-5,5])
view(2)
Pmax=max(max(P)) %presión máx y mín
Pmin=min(min(P))
}}
Matlab nos devuelve la siguiente imagen y los valores máximos y mínimos de la presión (Presión máx: 9.9959 Presión mínima: 6)
[[Archivo:Presionnn.jpg|800px|miniaturadeimagen|centro|Fig.PRESIÓN DEL FLUIDO]] 
Recordemos ahora la imagen del campo de velocidades (ahora sin ampliar) que habíamos obtenido anteriormente con el código matlab de el apartado 2.2: 

[[Archivo:velocidadd.jpg|800px|miniaturadeimagen|centro|Fig.CAMPO DE VELOCIDADES]] 

Si comparamos nuestra imagen de la presión del fluido (Fig.PRESIÓN DEL FLUIDO) con la del campo de velocidades (Fig.CAMPO DE VELOCIDADES) podemos observar que los puntos que cuentan con módulo de velocidad máxima (<math>\theta=0</math> y <math>\theta={\pi}</math>) son aquellos con presiones mínimas. Esto es un resultado coherente que verifica la ecuación de Bernoulli, a mayor velocidad, menores presiones.

== ºEcuación de Navier-Stokes ==

La ecuación de Navier-Stokes describe el movimiento de los fluidos a partir de su campo de velocidades y de presiones:

<math> d(\vec{u}·\nabla)\vec{u} + \nabla p = \mu \Delta \vec{u} </math> 

Los fluidos incompresibles tienen la propiedad de ocupar siempre el mismo volumen, de esta forma la viscosidad será nula: <math> \mu =0 </math> y la densidad será nula también: <math>d = 2</math>. Si sustituimos en la ecuación general, llegamos a la conclusión de que para un fluido incompresible:

<math> 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 </math> 

Si consideramos los campos <math> (\vec{u},p) </math> hallados anteriormente:

<math> (\vec{u}·\nabla)\vec{u} = u_{1} \frac{\partial u_{1}}{\partial \rho} \vec{e}_{\rho} + u_{1} \frac{\partial u_{1}}{\partial \theta} \vec{e}_{\theta}+ u_{2} \frac{\partial u_{2}}{\partial \rho} \vec{e}_{\rho}+ u_{2} \frac{\partial u_{2}}{\partial \theta} \vec{e}_{\theta} = </math> :

<math> (\vec{u}·\nabla)\vec{u}=(sen(\theta)(1-\frac{1}{\rho^2})[\frac{2}{\rho^3}sin(\theta)\vec{e}_{\rho}+(1-\frac{1}{\rho^2})cos(\theta)\vec{e}_{\theta}])+((1+\frac{1}{\rho^2})cos(\theta)[-\frac{2}{\rho^3}cos(\theta)\vec{e}_{\rho}-(1+\frac{1}{\rho^2})sin(\theta)\vec{e}_{\theta}])=</math>

<math> =(\vec{u}·\nabla)\vec{u}=[\frac{2}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)-\frac{2}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[(1-\frac{1}{\rho^2})^2sin(\theta)cos(\theta)-(1+\frac{1}{\rho^2})^2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

Lo aplicamos a nuestro caso con la densidad determinada:

<math> 2(\vec{u}·\nabla)\vec{u}=[\frac{4}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)- \frac{4}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[(1-\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)-(1+\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

De la misma forma:

<math> \nabla p = \frac{\partial p}{\partial \rho} \vec{e}_{\rho} + \frac{\partial p}{\partial \theta} \vec{e}_{\theta} = </math> :

<math> \nabla p = -[(2(1-\frac{1}{\rho^2})\frac{2}{\rho^3}sin^2(\theta))+(2(1+\frac{1}{\rho^2})(-\frac{2}{\rho^3})cos(\theta))]\vec{e}_{\rho} - [((1-\frac{1}{\rho^2})^2)2sin(\theta)cos(\theta)+((1+\frac{1}{\rho^2})^2)2(-sin(\theta))cos(\theta)]\vec{e}_{\theta}=</math>

<math> = [-\frac{4}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)+ \frac{4}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[-(1-\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)+(1+\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

Habiendo sustituido todos los valores en la ecuación de Navier-Stokes, coparando los resultados vemos directamente que se cumple <math> 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 </math> quedando así demostrado que se cumple la ecuación.

== º Trayectoria de la línea de corriente==
Si fuéramos una partícula del fluido con el que estamos trabajando, seguiríamos la trayectoria de una línea de corriente, es decir empezaríamos con su trayectoria original y a medida que se va acercando al obstáculo recorreríamos el contorno de este, y cuando lo hubiésemos pasado, volveríamos a la dirección inicial. En nuestro caso, como el obstáculo tiene un contorno circular el fluido empieza recto, luego se curva para evitar el obstáculo y finalmente vuelve a ir en la dirección inicial.
En relación a la variación de la velocidad y presión podemos comprobar con las gráficas, que en los puntos donde la velocidad es mínima las presiones son máximas. En los puntos que ya hemos calculado anteriormente donde la velocidad es 0, <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math>, sabemos que la presión es máxima.
Con esto podemos deducir que la presión y velocidad de estas partículas son inversamente proporcionales.

== º Paradoja D'Alembert==
:
<math>\vec{u}=sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}</math>

Para <math>\rho=1</math>:

<math>\vec u=2cos(\theta)\vec {e}_{\theta}</math>:

Según el Teorema de Kutta-Joukowski la fuerza que ejerce el fluido sobre el obstáculo, en este caso el círculo de radio 1, es proporcional a la circulación. Por lo tanto:

<math>\int_{0}^{2\pi} \vec u r'(\theta)\,d\theta =\int_{0}^{2\pi} 2cos(\theta)\vec{e}_{\theta}\vec{e}_{\theta}\,d\theta =\int_{0}^{2\pi} 2cos(\theta)\,d\theta=2[sin(0)-sin(2\pi)]=0 </math>:

El fluido no ejerce ninguna fuerza sobre el obstáculo, dando lugar a la paradoja D'Alembert, que concluye que la fuerza resultante sobre el cuerpo es cero, contradiciéndose con la observación. Esto demuestra también que la circulación del campo vectorial a lo largo de la circunferencia se anula.
Para poder ilustrar esta fuerza se necesitan tomar soluciones de modelos de fluidos con una viscosidad positiva.

== º Curvas de nivel de la presión==
Observando las curvas de nivel de la presión podemos entender mejor las conclusiones a las que hemos ido llegando a lo largo del trabajo. En esta gráfica podemos observar los puntos de presión máxima <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math> y los puntos de presión mínima <math>\theta=0</math> y <math>\theta=\pi</math>.
Además, cabe destacar que, en las zonas en las que las curvas de nivel están mas juntas, hay mayor variación de la presión.

[[Archivo:Doce_mnj.jpg|800px|miniaturadeimagen|centro|Líneas de presión]]

Código Matlab:
{{matlab| codigo=

contour(x,y,P,100) % Curvas de nivel de la presión
axis([-5,5,-5,5]) % Región a dibujar

}}
[[Categoría:Teoría de Campos]]
[[Categoría:TC21/22]]

Visualización de campos escalares y vectoriales en fluidos. (Grupo B6)

2021-12-11T15:21:18Z

Julia Cabeza: /* º Rotacional y Divergencia */

{{ TrabajoED | Visualización de campos escalares y vectoriales en fluidos. (Grupo B6) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC21/22|2021-22]] | María Magaldi Jurado, Nuria Martínez Ballester, Julia Cabeza Duque }}

En este trabajo vamos a visualizar campos escalares y vectoriales en fluidos. Más específicamente, el fluido a tratar es uno incompresible que se encuentra alrededor de un obstáculo de forma circular, formando así un anillo comprendido entre las circunferencias de radio 1 y 5, con el centro en el origen de coordenadas. Para facilitar el trabajo usaremos coordenadas polares (cilíndricas <math>(\rho,\theta)</math>). Para dibujar que el fluido va por el exterior de la superficie definida consideraremos que <math>(\rho,\theta) \in [-5,5]*[-5,5]</math>.
Además, tendremos que la velocidad de las partículas del fluido viene definida por el gradiente de la función potencial <math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math>

== º Mallado del fluido ==
El siguiente mallado representa los puntos interiores de la región ocupada por un fluido. Se encuentra en el exterior del círculo unidad. Como ya se ha introducido, nuestro mallado del anillo está comprendido entre los radios 1 y 5 y centro el origen.

[[Archivo:1,mallad0.jpg|800px|thumb|center|Mallado del fluido]]

El CÓDIGO matlab empleado es:
{{matlab|codigo=
ro=linspace(1,5,50); %creación de los intervalos
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt); %creación matrices
Mx=Mro.*cos(Mtt); %parametrización de x
My=Mro.*sin(Mtt); %parametrización de y
Mz=zeros(size(Mx)); %parametrización de z (nula al estar en un plano)
subplot(1,1,1);
mesh(Mx,My,Mz) %mallado
axis([-5,5,-5,5]) %superficie definida en [-5,5]*[-5,5]
view(2)
}}

== º Función potencial <math>\varphi</math> y velocidad de las partículas del fluido ==
=== º Función potencial ===
La velocidad de las partículas del fluido viene dada por el gradiente de la función potencial:
<math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math>
Las curvas de nivel de esta función potencial son las siguientes:
[[Archivo:2curvasdenivel.jpg|800px|thumb|center|FUNCIÓN POTENCIAL]]
El código matlab empleado es:
{{matlab|codigo=
ro=linspace(1,5,50);
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt);
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));
f=(Mro+1./Mro).*sin(Mtt); %función potencial
subplot(1,1,1)
contour(Mx,My,f,50); colorbar %curvas de nivel función potencial
}}

=== º Campo de velocidades (gradiente de función potencial) ===
Si dibujamos nuestro campo de velocidades <math>\vec u=\nabla\varphi</math> en matlab, la representación del mismo con flechas deberá ser ortogonal a las curvas de nivel antes dibujadas. Ya que, de forma geométrica, el gradiente es un vector normal a la curva de nivel de la función de la cual es gradiente en el punto que se esté estudiando.
Hemos pasado nuestra función potencial a cartesianas y trabajado posteriormente con sus derivadas. 
Así, con las igualdades: <math>\rho=\sqrt{x^2+y^2}\qquad\theta=arctg(\frac{y}{x})</math> 
Obtenemos la función potencial: <math>\varphi=y+\frac{y}{x^2+y^2}</math> 

El gradiente se define como:
<math>\vec u=\nabla\varphi=\frac{∂φ}{∂x}\vec i + \frac{∂φ}{∂j}\vec j </math>

Realizamos las derivadas parciales: 
<math>FX=\frac{\partial\varphi(x,y)}{\partial x}=\frac{-2xy}{(x^2+y^2)^2}\qquad FY=\frac{\partial\varphi(x,y)}{\partial y}=\frac{x^2-y^2}{(x^2+y^2)^2}</math> 
Con estas funciones, podremos dibujar en matlab el campo de velocidades con el siguiente código:
{{matlab|codigo=
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));
f=(Mro+1./Mro).*sin(Mtt); %función potencial
syms x y %nuevas variables en cartesianas para definir derivadas (campo de velocidades)
fx=inline('-2*(x.*y)./((x.^2+y.^2).^2)','x','y'); %derivadas que definen nuestro campos de velocidades
fy=inline('1+((x.^2-y.^2)./((x.^2+y.^2).^2))','x','y');
d1=fx(Mx,My);
d2=fy(Mx,My);
hold on
subplot(1,1,1)
contour(Mx,My,f,20); colorbar %curvas de nivel función potencial
quiver(Mx,My,d1,d2) %campo de velocidades fluido, ortogonal a curvas de nivel
axis([-5,5,-5,5])
hold off
}}
En nuestro dibujo saldrán las flechas, que representan nuestro campo de velocidades superpuestas con las curvas de nivel de la función potencial. Si ampliamos la imagen podremos observar que estos dos elementos son ortogonales:

[[Archivo: 2campoycurvas.jpg|800px|thumb|center|CAMPO DE VELOCIDADES]]

==º Ortogonalidad <math>\vec u</math> y <math>\vec n</math> ==
Como el obstáculo es un círculo situado en el plano xy, sabemos con certeza que entonces el vector normal será el correspondiente al eje z, es decir: <math>\vec n= \vec k</math>
Del apartado anterior conocemos que <math>\vec u=\nabla\varphi</math>.

<math>\vec{u}=sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}</math>

Conociendo entonces el vector normal y el vector gradiente, podemos calcular su producto escalar. Para comprobar si dos vectores son ortogonales, su producto escalar debe de ser nulo, es por ello que realizamos el producto escalar del gradiente y del vector normal.
<math>\nabla\varphi\cdot\vec n</math> = <math>(sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}})\cdot\vec ({e}_{z})</math> = 0

Como podemos observar, el producto escalar resulta 0, y queda comprobado así que el vector función potencial es perpendicular a su vector normal.

==º Calculo de <math>\vec u</math> con <math> \frac{1}{\rho}</math> despreciable ==
Al centrarnos en el fluido a una distancia muy lejana, podemos considerar <math> \frac{1}{\rho}</math> despreciable. En este caso sabremos que <math> u</math> cambia. La velocidad de las partículas desde los puntos lejanos del fluido viene definida por el gradiente de la función potencial <math>\varphi=\rho\sin(\theta)</math>

<math>\vec{u}=sin(\theta){\vec{e}_{\rho}}+\frac{1}{\rho}\rho\cos(\theta){\vec{e}_{\theta}}</math> 

En este caso:
<math>|\vec{u}| =\sqrt{cos(\theta)^2+sin(\theta)^2}=1</math>

== º Rotacional y Divergencia ==
:
<math>\nabla\times\vec u=\begin{vmatrix} \vec {e}_{\rho}&\rho\vec {e}_{\theta}&\vec {e}_{z} \\ \frac{\partial}{\partial{\rho}} & \frac{\partial}{\partial{\theta}} & \frac{\partial}{\partial{z}} \\ sin(\theta)[1-\frac{1}{\rho^2}] & \rho\cos(\theta)[\rho +\frac{1}{\rho}] & {0} \end{vmatrix}=[-cos(\theta)(1-\frac{1}{\rho^2})]\vec {e}_{z}+[cos(\theta)(1-\frac{1}{\rho^2})]\vec {e}_{z}=\vec {0}</math>
:
<math>\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}({\rho}sin(\theta)(1 -\frac{1}{\rho^2}))+\frac{\partial}{\partial{\theta}}(cos(\theta)(\rho +\frac{1}{\rho})+\frac{\partial}{\partial{z}}(\rho\cdot{0})]=\frac{1}{\rho}[sin(\theta)(1+\frac{1}{\rho^2})-sin(\theta)(1+\frac{1}{\rho^2})]=0 </math>

En este apartado comprobamos analíticamente que la divergencia y rotacional de <math> u</math> son nulas. Que la divergencia sea nula indica que en este caso el volumen del fluido no varía, es decir, se mantiene, ni se expande ni se contrae. A esta condición se le denomina incompresibilidad. Además, que el rotacional sea nulo demuestra que las partículas del fluido no giran.

== ºLíneas de corriente ==
Las líneas de corriente son aquellas que son tangentes a cada punto del campo <math>\vec u</math>. Para ello calcularemos el vector <math>\vec v=\vec k\times\vec u</math>

<math>\vec v=\frac{1}{\rho}\begin{vmatrix} \vec {e}_{\rho}&\rho\vec {e}_{\theta}&\vec {e}_{z} \\ {0} & {0} & {1} \\ sin(\theta)(1-\frac{1}{\rho^2}) & cos(\theta)(1+\frac{1}{\rho^2})\rho & {0} \end{vmatrix}= sin(\theta)(1-\frac{1}{\rho^2})\vec {e}_{\theta} - cos(\theta)(1+\frac{1}{\rho^2})\vec {e}_{\rho} =\vec {grad}\psi</math>:

<math>\psi=\int - cos(\theta)(1+\frac{1}{\rho^2}) \, d\rho \,\!=-cos(\theta)(\rho-\frac{1}{\rho}) + f(\theta) </math>:

<math>\frac{\partial\psi}{\partial\theta}=sin(\theta)(\rho-\frac{1}{\rho})+f'(\theta)=sin(\theta)(\rho-\frac{1}{\rho})\qquad f'(\theta)=0\qquad f(\theta)=cte</math>:

<math>\psi=-cos(\theta)(\rho-\frac{1}{\rho})</math>:

Del apartado anterior, podemos ver que el rotacional es nulo, y por lo tanto es irrotacional. De esta forma, las líneas de corriente serán tangentes al campo. 
Las curvas de nivel de la corriente son las siguientes: 
[[Archivo:curvascorriente.jpg|800px|miniaturadeimagen|centro|LÍNEAS DE CORRIENTE]] 
Si superponemos la imagen con el campo de velocidades antes calculado, observamos que efectivamente son tangentes: 
[[Archivo: curvasycampo.png|800px|miniaturadeimagen|centro|Representación de las líneas de corriente, tangentes al campo <math>\vec u</math>]]

El código Matlab empleado es el siguiente:

{{matlab|codigo=ro=linspace(1,5,300); %intervalo de ro (1,5)
tt=linspace(0,2*pi,300); %intervalo de teta (0,2pi)
[Mro,Mtt]=meshgrid(ro,tt);
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));

syms x y
fi=inline('-x+x./(x.^2+y.^2)','x','y'); %funciÃ³n que define la corriente de u
f2=fi(Mx,My);
fx=inline('-2*(x.*y)./((x.^2+y.^2).^2)','x','y'); %derivadas de funcion potencial que definen campo de velocidades del fluido
fy=inline('1+((x.^2-y.^2)./((x.^2+y.^2).^2))','x','y');
d1=fx(Mx,My);
d2=fy(Mx,My);
hold on
subplot(1,1,1)
contour(Mx,My,f2,100) %curvas de nivel de la corriente
quiver(Mx,My,d1,d2) %campo de velocidades fluido, tangente a lineas de corriente
hold off
}}

== º Puntos frontera==
Al tener definida la superficie del anillo con los radios 1 y 5, buscamos donde están situados los puntos de frontera. Con la función potencial inicial <math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math> y con
<math> u</math> = <math>(sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}) </math>.

Calculamos el modulo del vector en función de <math>\theta</math> para

<math>\rho=1\qquad </math> <math> u</math> = <math> 2cos(\theta){\vec{e}_{\theta}} </math>

<math>|\vec u|=\sqrt{4cos(\theta)^2} = 2cos(\theta) = 2|cos(\theta)|</math> 

Como podemos observar, el módulo de la velocidad es máximo cuando cos<math>\theta</math> alcanza sus valores máximos (en 0 y en <math>\pi</math>, donde alcanza en valor absoluto el valor de 1). 

Así, los puntos de velocidad máxima son <math>\theta=0</math> y <math>\theta={pi}</math>. 

Los puntos de velocidad mínima o puntos de remanso son <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math>. 

== º Presión del fluido ==

En dinámica de fluidos, el principio de Bernoulli, también denominado ecuación de Bernoulli, describe el comportamiento de un fluido moviéndose a lo largo de una línea de corriente. Relacionando las siguientes variables: 

<math>p =</math>presión estática a la que está sometido el fluido debida a las moléculas que lo rodean. 
<math>d =</math>densidad del fluido. 

<math>\vec { u } =</math>velocidad de flujo del fluido. :

<math>\frac { 1 }{ 2 } d { \left| \vec { u } \right| }^{ 2 }+ { p } =cte</math> 
Para que esta expresión pueda aplicarse, el fluido debe ser incompresible, no viscoso y fluir en régimen estacionario (la velocidad en un punto determinado no varía con el tiempo). Tomaremos nuestro fluido como tal. 
Suponiendo que nuestro fluido efectivamente cumple la ecuación de Bernouilli, que posee una densidad constante de d=2 y tomando la cte de la ecuación como cte=10, calcularemos la presión de nuestro fluido. 
Sustituyendo los valores mencionados tenemos que: 
<math> p = 10- |\vec{u}|^2 </math> 
Calculamos el módulo al cuadrado de nuestro campo de velocidades: 
<math> |\vec{u}|^2 = (1-\frac{1}{ρ^2})^2 sen^2(θ) + (1+\frac{1}{ρ^2})^2 cos^2(θ) </math> 
Así la función que define la presión del fluido es: 
<math> p = 10- (1-\frac{1}{ρ^2})^2 sen^2(θ) + (1+\frac{1}{ρ^2})^2 cos^2(θ) </math> 

Para la representación de la presión empleamos el siguiente código:
{{matlab| codigo=
ro=linspace(1,5,50);
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt);
x=Mro.*cos(Mtt);
y=Mro.*sin(Mtt);
P=10-(((1-1./Mro.^2).^2).*(sin(Mtt).^2)+((1+1./Mro.^2).^2).*(cos(Mtt).^2)); %función presión
surf(x,y,P); colorbar;
axis([-5,5,-5,5])
view(2)
Pmax=max(max(P)) %presión máx y mín
Pmin=min(min(P))
}}
Matlab nos devuelve la siguiente imagen y los valores máximos y mínimos de la presión (Presión máx: 9.9959 Presión mínima: 6)
[[Archivo:Presionnn.jpg|800px|miniaturadeimagen|centro|Fig.PRESIÓN DEL FLUIDO]] 
Recordemos ahora la imagen del campo de velocidades (ahora sin ampliar) que habíamos obtenido anteriormente con el código matlab de el apartado 2.2: 

[[Archivo:velocidadd.jpg|800px|miniaturadeimagen|centro|Fig.CAMPO DE VELOCIDADES]] 

Si comparamos nuestra imagen de la presión del fluido (Fig.PRESIÓN DEL FLUIDO) con la del campo de velocidades (Fig.CAMPO DE VELOCIDADES) podemos observar que los puntos que cuentan con módulo de velocidad máxima (<math>\theta=0</math> y <math>\theta={\pi}</math>) son aquellos con presiones mínimas. Esto es un resultado coherente que verifica la ecuación de Bernoulli, a mayor velocidad, menores presiones.

== ºEcuación de Navier-Stokes ==

La ecuación de Navier-Stokes describe el movimiento de los fluidos a partir de su campo de velocidades y de presiones:

<math> d(\vec{u}·\nabla)\vec{u} + \nabla p = \mu \Delta \vec{u} </math> 

Los fluidos incompresibles tienen la propiedad de ocupar siempre el mismo volumen, de esta forma la viscosidad será nula: <math> \mu =0 </math> y la densidad será nula también: <math>d = 2</math>. Si sustituimos en la ecuación general, llegamos a la conclusión de que para un fluido incompresible:

<math> 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 </math> 

Si consideramos los campos <math> (\vec{u},p) </math> hallados anteriormente:

<math> (\vec{u}·\nabla)\vec{u} = u_{1} \frac{\partial u_{1}}{\partial \rho} \vec{e}_{\rho} + u_{1} \frac{\partial u_{1}}{\partial \theta} \vec{e}_{\theta}+ u_{2} \frac{\partial u_{2}}{\partial \rho} \vec{e}_{\rho}+ u_{2} \frac{\partial u_{2}}{\partial \theta} \vec{e}_{\theta} = </math> :

<math> (\vec{u}·\nabla)\vec{u}=(sen(\theta)(1-\frac{1}{\rho^2})[\frac{2}{\rho^3}sin(\theta)\vec{e}_{\rho}+(1-\frac{1}{\rho^2})cos(\theta)\vec{e}_{\theta}])+((1+\frac{1}{\rho^2})cos(\theta)[-\frac{2}{\rho^3}cos(\theta)\vec{e}_{\rho}-(1+\frac{1}{\rho^2})sin(\theta)\vec{e}_{\theta}])=</math>

<math> =(\vec{u}·\nabla)\vec{u}=[\frac{2}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)-\frac{2}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[(1-\frac{1}{\rho^2})^2sin(\theta)cos(\theta)-(1+\frac{1}{\rho^2})^2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

Lo aplicamos a nuestro caso con la densidad determinada:

<math> 2(\vec{u}·\nabla)\vec{u}=[\frac{4}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)- \frac{4}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[(1-\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)-(1+\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

De la misma forma:

<math> \nabla p = \frac{\partial p}{\partial \rho} \vec{e}_{\rho} + \frac{\partial p}{\partial \theta} \vec{e}_{\theta} = </math> :

<math> \nabla p = -[(2(1-\frac{1}{\rho^2})\frac{2}{\rho^3}sin^2(\theta))+(2(1+\frac{1}{\rho^2})(-\frac{2}{\rho^3})cos(\theta))]\vec{e}_{\rho} - [((1-\frac{1}{\rho^2})^2)2sin(\theta)cos(\theta)+((1+\frac{1}{\rho^2})^2)2(-sin(\theta))cos(\theta)]\vec{e}_{\theta}=</math>

<math> = [-\frac{4}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)+ \frac{4}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[-(1-\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)+(1+\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

Habiendo sustituido todos los valores en la ecuación de Navier-Stokes, coparando los resultados vemos directamente que se cumple <math> 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 </math> quedando así demostrado que se cumple la ecuación.

== º Trayectoria de la línea de corriente==
Si fuéramos una partícula del fluido con el que estamos trabajando, seguiríamos la trayectoria de una línea de corriente, es decir empezaríamos con su trayectoria original y a medida que se va acercando al obstáculo recorreríamos el contorno de este, y cuando lo hubiésemos pasado terminado volveríamos a la dirección inicial. En nuestro caso, como el obstáculo tiene un contorno circular el fluido empieza recto, luego se curva para evitar el obstáculo y finalmente vuelve a ir en la dirección inicial.
En relación a la variación de la velocidad y presión podemos comprobar con las gráficas, que en los puntos donde la velocidad es mínima las presiones son máximas. En los puntos que ya hemos calculado anteriormente donde la velocidad es 0, <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math>, sabemos que la presión es máxima.
Con esto podemos deducir que la presión y velocidad de estas partículas son inversamente proporcionales.

== º Paradoja D'Alembert==
:
<math>\vec{u}=sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}</math>

Para <math>\rho=1</math>:

<math>\vec u=2cos(\theta)\vec {e}_{\theta}</math>:

Según el Teorema de Kutta-Joukowski la fuerza que ejerce el fluido sobre el obstáculo, en este caso el círculo de radio 1, es proporcional a la circulación. Por lo tanto:

<math>\int_{0}^{2\pi} \vec u r'(\theta)\,d\theta =\int_{0}^{2\pi} 2cos(\theta)\vec{e}_{\theta}\vec{e}_{\theta}\,d\theta =\int_{0}^{2\pi} 2cos(\theta)\,d\theta=2[sin(0)-sin(2\pi)]=0 </math>:

El fluido no ejerce ninguna fuerza sobre el obstáculo, dando lugar a la paradoja D'Alembert, que concluye que la fuerza resultante sobre el cuerpo es cero, contradiciéndose con la observación. Esto demuestra también que la circulación del campo vectorial a lo largo de la circunferencia se anula.
Para poder ilustrar esta fuerza se necesitan tomar soluciones de modelos de fluidos con una viscosidad positiva.

== º Curvas de nivel de la presión==
Observando las curvas de nivel de la presión podemos entender mejor las conclusiones a las que hemos ido llegando a lo largo del trabajo. En esta gráfica podemos observar los puntos de presión máxima <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math> y los puntos de presión mínima <math>\theta=0</math> y <math>\theta=\pi</math>.
Además, cabe destacar que, en las zonas en las que las curvas de nivel están mas juntas, hay mayor variación de la presión.

[[Archivo:Doce_mnj.jpg|800px|miniaturadeimagen|centro|Líneas de presión]]

Código Matlab:
{{matlab| codigo=

contour(x,y,P,100) % Curvas de nivel de la presión
axis([-5,5,-5,5]) % Región a dibujar

}}
[[Categoría:Teoría de Campos]]
[[Categoría:TC21/22]]

Visualización de campos escalares y vectoriales en fluidos. (Grupo B6)

2021-12-11T15:14:31Z

Julia Cabeza: /* º Ortogonalidad \nabla\vec u y \vec n */

{{ TrabajoED | Visualización de campos escalares y vectoriales en fluidos. (Grupo B6) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC21/22|2021-22]] | María Magaldi Jurado, Nuria Martínez Ballester, Julia Cabeza Duque }}

En este trabajo vamos a visualizar campos escalares y vectoriales en fluidos. Más específicamente, el fluido a tratar es uno incompresible que se encuentra alrededor de un obstáculo de forma circular, formando así un anillo comprendido entre las circunferencias de radio 1 y 5, con el centro en el origen de coordenadas. Para facilitar el trabajo usaremos coordenadas polares (cilíndricas <math>(\rho,\theta)</math>). Para dibujar que el fluido va por el exterior de la superficie definida consideraremos que <math>(\rho,\theta) \in [-5,5]*[-5,5]</math>.
Además, tendremos que la velocidad de las partículas del fluido viene definida por el gradiente de la función potencial <math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math>

== º Mallado del fluido ==
El siguiente mallado representa los puntos interiores de la región ocupada por un fluido. Se encuentra en el exterior del círculo unidad. Como ya se ha introducido, nuestro mallado del anillo está comprendido entre los radios 1 y 5 y centro el origen.

[[Archivo:1,mallad0.jpg|800px|thumb|center|Mallado del fluido]]

El CÓDIGO matlab empleado es:
{{matlab|codigo=
ro=linspace(1,5,50); %creación de los intervalos
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt); %creación matrices
Mx=Mro.*cos(Mtt); %parametrización de x
My=Mro.*sin(Mtt); %parametrización de y
Mz=zeros(size(Mx)); %parametrización de z (nula al estar en un plano)
subplot(1,1,1);
mesh(Mx,My,Mz) %mallado
axis([-5,5,-5,5]) %superficie definida en [-5,5]*[-5,5]
view(2)
}}

== º Función potencial <math>\varphi</math> y velocidad de las partículas del fluido ==
=== º Función potencial ===
La velocidad de las partículas del fluido viene dada por el gradiente de la función potencial:
<math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math>
Las curvas de nivel de esta función potencial son las siguientes:
[[Archivo:2curvasdenivel.jpg|800px|thumb|center|FUNCIÓN POTENCIAL]]
El código matlab empleado es:
{{matlab|codigo=
ro=linspace(1,5,50);
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt);
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));
f=(Mro+1./Mro).*sin(Mtt); %función potencial
subplot(1,1,1)
contour(Mx,My,f,50); colorbar %curvas de nivel función potencial
}}

=== º Campo de velocidades (gradiente de función potencial) ===
Si dibujamos nuestro campo de velocidades <math>\vec u=\nabla\varphi</math> en matlab, la representación del mismo con flechas deberá ser ortogonal a las curvas de nivel antes dibujadas. Ya que, de forma geométrica, el gradiente es un vector normal a la curva de nivel de la función de la cual es gradiente en el punto que se esté estudiando.
Hemos pasado nuestra función potencial a cartesianas y trabajado posteriormente con sus derivadas. 
Así, con las igualdades: <math>\rho=\sqrt{x^2+y^2}\qquad\theta=arctg(\frac{y}{x})</math> 
Obtenemos la función potencial: <math>\varphi=y+\frac{y}{x^2+y^2}</math> 

El gradiente se define como:
<math>\vec u=\nabla\varphi=\frac{∂φ}{∂x}\vec i + \frac{∂φ}{∂j}\vec j </math>

Realizamos las derivadas parciales: 
<math>FX=\frac{\partial\varphi(x,y)}{\partial x}=\frac{-2xy}{(x^2+y^2)^2}\qquad FY=\frac{\partial\varphi(x,y)}{\partial y}=\frac{x^2-y^2}{(x^2+y^2)^2}</math> 
Con estas funciones, podremos dibujar en matlab el campo de velocidades con el siguiente código:
{{matlab|codigo=
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));
f=(Mro+1./Mro).*sin(Mtt); %función potencial
syms x y %nuevas variables en cartesianas para definir derivadas (campo de velocidades)
fx=inline('-2*(x.*y)./((x.^2+y.^2).^2)','x','y'); %derivadas que definen nuestro campos de velocidades
fy=inline('1+((x.^2-y.^2)./((x.^2+y.^2).^2))','x','y');
d1=fx(Mx,My);
d2=fy(Mx,My);
hold on
subplot(1,1,1)
contour(Mx,My,f,20); colorbar %curvas de nivel función potencial
quiver(Mx,My,d1,d2) %campo de velocidades fluido, ortogonal a curvas de nivel
axis([-5,5,-5,5])
hold off
}}
En nuestro dibujo saldrán las flechas, que representan nuestro campo de velocidades superpuestas con las curvas de nivel de la función potencial. Si ampliamos la imagen podremos observar que estos dos elementos son ortogonales:

[[Archivo: 2campoycurvas.jpg|800px|thumb|center|CAMPO DE VELOCIDADES]]

==º Ortogonalidad <math>\vec u</math> y <math>\vec n</math> ==
Como el obstáculo es un círculo situado en el plano xy, sabemos con certeza que entonces el vector normal será el correspondiente al eje z, es decir: <math>\vec n= \vec k</math>
Del apartado anterior conocemos que <math>\vec u=\nabla\varphi</math>.

<math>\vec{u}=sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}</math>

Conociendo entonces el vector normal y el vector gradiente, podemos calcular su producto escalar. Para comprobar si dos vectores son ortogonales, su producto escalar debe de ser nulo, es por ello que realizamos el producto escalar del gradiente y del vector normal.
<math>\nabla\varphi\cdot\vec n</math> = <math>(sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}})\cdot\vec ({e}_{z})</math> = 0

Como podemos observar, el producto escalar resulta 0, y queda comprobado así que el vector función potencial es perpendicular a su vector normal.

==º Calculo de <math>\vec u</math> con <math> \frac{1}{\rho}</math> despreciable ==
Al centrarnos en el fluido a una distancia muy lejana, podemos considerar <math> \frac{1}{\rho}</math> despreciable. En este caso sabremos que <math> u</math> cambia. La velocidad de las partículas desde los puntos lejanos del fluido viene definida por el gradiente de la función potencial <math>\varphi=\rho\sin(\theta)</math>

<math>\vec{u}=sin(\theta){\vec{e}_{\rho}}+\frac{1}{\rho}\rho\cos(\theta){\vec{e}_{\theta}}</math> 

En este caso:
<math>|\vec{u}| =\sqrt{cos(\theta)^2+sin(\theta)^2}=1</math>

== º Rotacional y Divergencia ==
:
<math>\nabla\times\vec u=\begin{vmatrix} \vec {e}_{\rho}&\vec {e}_{\theta}&\vec {e}_{z} \\ \frac{\partial}{\partial{\rho}} & \frac{\partial}{\partial{\theta}} & \frac{\partial}{\partial{z}} \\ sin(\theta)[1-\frac{1}{\rho^2}] & cos(\theta)[\rho +\frac{1}{\rho}] & {0} \end{vmatrix}=[-cos(\theta)(1-\frac{1}{\rho^2})]\vec {e}_{z}+[cos(\theta)(1-\frac{1}{\rho^2})]\vec {e}_{z}=\vec {0}</math>
:
<math>\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}({\rho}sin(\theta)(1 -\frac{1}{\rho^2}))+\frac{\partial}{\partial{\theta}}(cos(\theta)(\rho +\frac{1}{\rho})+\frac{\partial}{\partial{z}}(\rho\cdot{0})]=\frac{1}{\rho}[sin(\theta)(1+\frac{1}{\rho^2})-sin(\theta)(1+\frac{1}{\rho^2})]=0 </math>

En este apartado comprobamos analíticamente que la divergencia y rotacional de <math> u</math> son nulas. Que la divergencia sea nula indica que en este caso el volumen del fluido no varía, es decir, se mantiene, ni se expande ni se contrae. A esta condición se le denomina incompresibilidad. Además, que el rotacional sea nulo demuestra que las partículas del fluido no giran.

== ºLíneas de corriente ==
Las líneas de corriente son aquellas que son tangentes a cada punto del campo <math>\vec u</math>. Para ello calcularemos el vector <math>\vec v=\vec k\times\vec u</math>

<math>\vec v=\frac{1}{\rho}\begin{vmatrix} \vec {e}_{\rho}&\rho\vec {e}_{\theta}&\vec {e}_{z} \\ {0} & {0} & {1} \\ sin(\theta)(1-\frac{1}{\rho^2}) & cos(\theta)(1+\frac{1}{\rho^2})\rho & {0} \end{vmatrix}= sin(\theta)(1-\frac{1}{\rho^2})\vec {e}_{\theta} - cos(\theta)(1+\frac{1}{\rho^2})\vec {e}_{\rho} =\vec {grad}\psi</math>:

<math>\psi=\int - cos(\theta)(1+\frac{1}{\rho^2}) \, d\rho \,\!=-cos(\theta)(\rho-\frac{1}{\rho}) + f(\theta) </math>:

<math>\frac{\partial\psi}{\partial\theta}=sin(\theta)(\rho-\frac{1}{\rho})+f'(\theta)=sin(\theta)(\rho-\frac{1}{\rho})\qquad f'(\theta)=0\qquad f(\theta)=cte</math>:

<math>\psi=-cos(\theta)(\rho-\frac{1}{\rho})</math>:

Del apartado anterior, podemos ver que el rotacional es nulo, y por lo tanto es irrotacional. De esta forma, las líneas de corriente serán tangentes al campo. 
Las curvas de nivel de la corriente son las siguientes: 
[[Archivo:curvascorriente.jpg|800px|miniaturadeimagen|centro|LÍNEAS DE CORRIENTE]] 
Si superponemos la imagen con el campo de velocidades antes calculado, observamos que efectivamente son tangentes: 
[[Archivo: curvasycampo.png|800px|miniaturadeimagen|centro|Representación de las líneas de corriente, tangentes al campo <math>\vec u</math>]]

El código Matlab empleado es el siguiente:

{{matlab|codigo=ro=linspace(1,5,300); %intervalo de ro (1,5)
tt=linspace(0,2*pi,300); %intervalo de teta (0,2pi)
[Mro,Mtt]=meshgrid(ro,tt);
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));

syms x y
fi=inline('-x+x./(x.^2+y.^2)','x','y'); %funciÃ³n que define la corriente de u
f2=fi(Mx,My);
fx=inline('-2*(x.*y)./((x.^2+y.^2).^2)','x','y'); %derivadas de funcion potencial que definen campo de velocidades del fluido
fy=inline('1+((x.^2-y.^2)./((x.^2+y.^2).^2))','x','y');
d1=fx(Mx,My);
d2=fy(Mx,My);
hold on
subplot(1,1,1)
contour(Mx,My,f2,100) %curvas de nivel de la corriente
quiver(Mx,My,d1,d2) %campo de velocidades fluido, tangente a lineas de corriente
hold off
}}

== º Puntos frontera==
Al tener definida la superficie del anillo con los radios 1 y 5, buscamos donde están situados los puntos de frontera. Con la función potencial inicial <math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math> y con
<math> u</math> = <math>(sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}) </math>.

Calculamos el modulo del vector en función de <math>\theta</math> para

<math>\rho=1\qquad </math> <math> u</math> = <math> 2cos(\theta){\vec{e}_{\theta}} </math>

<math>|\vec u|=\sqrt{4cos(\theta)^2} = 2cos(\theta) = 2|cos(\theta)|</math> 

Como podemos observar, el módulo de la velocidad es máximo cuando cos<math>\theta</math> alcanza sus valores máximos (en 0 y en <math>\pi</math>, donde alcanza en valor absoluto el valor de 1). 

Así, los puntos de velocidad máxima son <math>\theta=0</math> y <math>\theta={pi}</math>. 

Los puntos de velocidad mínima o puntos de remanso son <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math>. 

== º Presión del fluido ==

En dinámica de fluidos, el principio de Bernoulli, también denominado ecuación de Bernoulli, describe el comportamiento de un fluido moviéndose a lo largo de una línea de corriente. Relacionando las siguientes variables: 

<math>p =</math>presión estática a la que está sometido el fluido debida a las moléculas que lo rodean. 
<math>d =</math>densidad del fluido. 

<math>\vec { u } =</math>velocidad de flujo del fluido. :

<math>\frac { 1 }{ 2 } d { \left| \vec { u } \right| }^{ 2 }+ { p } =cte</math> 
Para que esta expresión pueda aplicarse, el fluido debe ser incompresible, no viscoso y fluir en régimen estacionario (la velocidad en un punto determinado no varía con el tiempo). Tomaremos nuestro fluido como tal. 
Suponiendo que nuestro fluido efectivamente cumple la ecuación de Bernouilli, que posee una densidad constante de d=2 y tomando la cte de la ecuación como cte=10, calcularemos la presión de nuestro fluido. 
Sustituyendo los valores mencionados tenemos que: 
<math> p = 10- |\vec{u}|^2 </math> 
Calculamos el módulo al cuadrado de nuestro campo de velocidades: 
<math> |\vec{u}|^2 = (1-\frac{1}{ρ^2})^2 sen^2(θ) + (1+\frac{1}{ρ^2})^2 cos^2(θ) </math> 
Así la función que define la presión del fluido es: 
<math> p = 10- (1-\frac{1}{ρ^2})^2 sen^2(θ) + (1+\frac{1}{ρ^2})^2 cos^2(θ) </math> 

Para la representación de la presión empleamos el siguiente código:
{{matlab| codigo=
ro=linspace(1,5,50);
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt);
x=Mro.*cos(Mtt);
y=Mro.*sin(Mtt);
P=10-(((1-1./Mro.^2).^2).*(sin(Mtt).^2)+((1+1./Mro.^2).^2).*(cos(Mtt).^2)); %función presión
surf(x,y,P); colorbar;
axis([-5,5,-5,5])
view(2)
Pmax=max(max(P)) %presión máx y mín
Pmin=min(min(P))
}}
Matlab nos devuelve la siguiente imagen y los valores máximos y mínimos de la presión (Presión máx: 9.9959 Presión mínima: 6)
[[Archivo:Presionnn.jpg|800px|miniaturadeimagen|centro|Fig.PRESIÓN DEL FLUIDO]] 
Recordemos ahora la imagen del campo de velocidades (ahora sin ampliar) que habíamos obtenido anteriormente con el código matlab de el apartado 2.2: 

[[Archivo:velocidadd.jpg|800px|miniaturadeimagen|centro|Fig.CAMPO DE VELOCIDADES]] 

Si comparamos nuestra imagen de la presión del fluido (Fig.PRESIÓN DEL FLUIDO) con la del campo de velocidades (Fig.CAMPO DE VELOCIDADES) podemos observar que los puntos que cuentan con módulo de velocidad máxima (<math>\theta=0</math> y <math>\theta={\pi}</math>) son aquellos con presiones mínimas. Esto es un resultado coherente que verifica la ecuación de Bernoulli, a mayor velocidad, menores presiones.

== ºEcuación de Navier-Stokes ==

La ecuación de Navier-Stokes describe el movimiento de los fluidos a partir de su campo de velocidades y de presiones:

<math> d(\vec{u}·\nabla)\vec{u} + \nabla p = \mu \Delta \vec{u} </math> 

Los fluidos incompresibles tienen la propiedad de ocupar siempre el mismo volumen, de esta forma la viscosidad será nula: <math> \mu =0 </math> y la densidad será nula también: <math>d = 2</math>. Si sustituimos en la ecuación general, llegamos a la conclusión de que para un fluido incompresible:

<math> 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 </math> 

Si consideramos los campos <math> (\vec{u},p) </math> hallados anteriormente:

<math> (\vec{u}·\nabla)\vec{u} = u_{1} \frac{\partial u_{1}}{\partial \rho} \vec{e}_{\rho} + u_{1} \frac{\partial u_{1}}{\partial \theta} \vec{e}_{\theta}+ u_{2} \frac{\partial u_{2}}{\partial \rho} \vec{e}_{\rho}+ u_{2} \frac{\partial u_{2}}{\partial \theta} \vec{e}_{\theta} = </math> :

<math> (\vec{u}·\nabla)\vec{u}=(sen(\theta)(1-\frac{1}{\rho^2})[\frac{2}{\rho^3}sin(\theta)\vec{e}_{\rho}+(1-\frac{1}{\rho^2})cos(\theta)\vec{e}_{\theta}])+((1+\frac{1}{\rho^2})cos(\theta)[-\frac{2}{\rho^3}cos(\theta)\vec{e}_{\rho}-(1+\frac{1}{\rho^2})sin(\theta)\vec{e}_{\theta}])=</math>

<math> =(\vec{u}·\nabla)\vec{u}=[\frac{2}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)-\frac{2}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[(1-\frac{1}{\rho^2})^2sin(\theta)cos(\theta)-(1+\frac{1}{\rho^2})^2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

Lo aplicamos a nuestro caso con la densidad determinada:

<math> 2(\vec{u}·\nabla)\vec{u}=[\frac{4}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)- \frac{4}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[(1-\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)-(1+\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

De la misma forma:

<math> \nabla p = \frac{\partial p}{\partial \rho} \vec{e}_{\rho} + \frac{\partial p}{\partial \theta} \vec{e}_{\theta} = </math> :

<math> \nabla p = -[(2(1-\frac{1}{\rho^2})\frac{2}{\rho^3}sin^2(\theta))+(2(1+\frac{1}{\rho^2})(-\frac{2}{\rho^3})cos(\theta))]\vec{e}_{\rho} - [((1-\frac{1}{\rho^2})^2)2sin(\theta)cos(\theta)+((1+\frac{1}{\rho^2})^2)2(-sin(\theta))cos(\theta)]\vec{e}_{\theta}=</math>

<math> = [-\frac{4}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)+ \frac{4}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[-(1-\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)+(1+\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

Habiendo sustituido todos los valores en la ecuación de Navier-Stokes, coparando los resultados vemos directamente que se cumple <math> 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 </math> quedando así demostrado que se cumple la ecuación.

== º Trayectoria de la línea de corriente==
Si fuéramos una partícula del fluido con el que estamos trabajando, seguiríamos la trayectoria de una línea de corriente, es decir empezaríamos con su trayectoria original y a medida que se va acercando al obstáculo recorreríamos el contorno de este, y cuando lo hubiésemos pasado terminado volveríamos a la dirección inicial. En nuestro caso, como el obstáculo tiene un contorno circular el fluido empieza recto, luego se curva para evitar el obstáculo y finalmente vuelve a ir en la dirección inicial.
En relación a la variación de la velocidad y presión podemos comprobar con las gráficas, que en los puntos donde la velocidad es mínima las presiones son máximas. En los puntos que ya hemos calculado anteriormente donde la velocidad es 0, <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math>, sabemos que la presión es máxima.
Con esto podemos deducir que la presión y velocidad de estas partículas son inversamente proporcionales.

== º Paradoja D'Alembert==
:
<math>\vec{u}=sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}</math>

Para <math>\rho=1</math>:

<math>\vec u=2cos(\theta)\vec {e}_{\theta}</math>:

Según el Teorema de Kutta-Joukowski la fuerza que ejerce el fluido sobre el obstáculo, en este caso el círculo de radio 1, es proporcional a la circulación. Por lo tanto:

<math>\int_{0}^{2\pi} \vec u r'(\theta)\,d\theta =\int_{0}^{2\pi} 2cos(\theta)\vec{e}_{\theta}\vec{e}_{\theta}\,d\theta =\int_{0}^{2\pi} 2cos(\theta)\,d\theta=2[sin(0)-sin(2\pi)]=0 </math>:

El fluido no ejerce ninguna fuerza sobre el obstáculo, dando lugar a la paradoja D'Alembert, que concluye que la fuerza resultante sobre el cuerpo es cero, contradiciéndose con la observación. Esto demuestra también que la circulación del campo vectorial a lo largo de la circunferencia se anula.
Para poder ilustrar esta fuerza se necesitan tomar soluciones de modelos de fluidos con una viscosidad positiva.

== º Curvas de nivel de la presión==
Observando las curvas de nivel de la presión podemos entender mejor las conclusiones a las que hemos ido llegando a lo largo del trabajo. En esta gráfica podemos observar los puntos de presión máxima <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math> y los puntos de presión mínima <math>\theta=0</math> y <math>\theta=\pi</math>.
Además, cabe destacar que, en las zonas en las que las curvas de nivel están mas juntas, hay mayor variación de la presión.

[[Archivo:Doce_mnj.jpg|800px|miniaturadeimagen|centro|Líneas de presión]]

Código Matlab:
{{matlab| codigo=

contour(x,y,P,100) % Curvas de nivel de la presión
axis([-5,5,-5,5]) % Región a dibujar

}}
[[Categoría:Teoría de Campos]]
[[Categoría:TC21/22]]

Visualización de campos escalares y vectoriales en fluidos. (Grupo B6)

2021-12-10T11:35:20Z

Julia Cabeza: /* ºLíneas de corriente */

{{ TrabajoED | Visualización de campos escalares y vectoriales en fluidos. (Grupo B6) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC21/22|2021-22]] | María Magaldi Jurado, Nuria Martínez Ballester, Julia Cabeza Duque }}

En este trabajo vamos a visualizar campos escalares y vectoriales en fluidos. Más específicamente, el fluido a tratar es uno incompresible que se encuentra alrededor de un obstáculo de forma circular, formando así un anillo comprendido entre las circunferencias de radio 1 y 5, con el centro en el origen de coordenadas. Para facilitar el trabajo usaremos coordenadas polares (cilíndricas <math>(\rho,\theta)</math>). Para dibujar que el fluido va por el exterior de la superficie definida consideraremos que <math>(\rho,\theta) \in [-5,5]*[-5,5]</math>.
Además, tendremos que la velocidad de las partículas del fluido viene definida por el gradiente de la función potencial <math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math>

== º Mallado del fluido ==
El siguiente mallado representa los puntos interiores de la región ocupada por un fluido. Se encuentra en el exterior del círculo unidad. Como ya se ha introducido, nuestro mallado del anillo está comprendido entre los radios 1 y 5 y centro el origen.

[[Archivo:1,mallad0.jpg|800px|thumb|center|Mallado del fluido]]

El CÓDIGO matlab empleado es:
{{matlab|codigo=
ro=linspace(1,5,50); %creación de los intervalos
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt); %creación matrices
Mx=Mro.*cos(Mtt); %parametrización de x
My=Mro.*sin(Mtt); %parametrización de y
Mz=zeros(size(Mx)); %parametrización de z (nula al estar en un plano)
subplot(1,1,1);
mesh(Mx,My,Mz) %mallado
axis([-5,5,-5,5]) %superficie definida en [-5,5]*[-5,5]
view(2)
}}

== º Función potencial <math>\varphi</math> y velocidad de las partículas del fluido ==
=== º Función potencial ===
La velocidad de las partículas del fluido viene dada por el gradiente de la función potencial:
<math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math>
Las curvas de nivel de esta función potencial son las siguientes:
[[Archivo:2curvasdenivel.jpg|800px|thumb|center|FUNCIÓN POTENCIAL]]
El código matlab empleado es:
{{matlab|codigo=
ro=linspace(1,5,50);
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt);
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));
f=(Mro+1./Mro).*sin(Mtt); %función potencial
subplot(1,1,1)
contour(Mx,My,f,50); colorbar %curvas de nivel función potencial
}}

=== º Campo de velocidades (gradiente de función potencial) ===
Si dibujamos nuestro campo de velocidades <math>\vec u=\nabla\varphi</math> en matlab, la representación del mismo con flechas deberá ser ortogonal a las curvas de nivel antes dibujadas. Ya que, de forma geométrica, el gradiente es un vector normal a la curva de nivel de la función de la cual es gradiente en el punto que se esté estudiando.
Hemos pasado nuestra función potencial a cartesianas y trabajado posteriormente con sus derivadas. 
Así, con las igualdades: <math>\rho=\sqrt{x^2+y^2}\qquad\theta=arctg(\frac{y}{x})</math> 
Obtenemos la función potencial: <math>\varphi=y+\frac{y}{x^2+y^2}</math> 

El gradiente se define como:
<math>\vec u=\nabla\varphi=\frac{∂φ}{∂x}\vec i + \frac{∂φ}{∂j}\vec j </math>

Realizamos las derivadas parciales: 
<math>FX=\frac{\partial\varphi(x,y)}{\partial x}=\frac{-2xy}{(x^2+y^2)^2}\qquad FY=\frac{\partial\varphi(x,y)}{\partial y}=\frac{x^2-y^2}{(x^2+y^2)^2}</math> 
Con estas funciones, podremos dibujar en matlab el campo de velocidades con el siguiente código:
{{matlab|codigo=
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));
f=(Mro+1./Mro).*sin(Mtt); %función potencial
syms x y %nuevas variables en cartesianas para definir derivadas (campo de velocidades)
fx=inline('-2*(x.*y)./((x.^2+y.^2).^2)','x','y'); %derivadas que definen nuestro campos de velocidades
fy=inline('1+((x.^2-y.^2)./((x.^2+y.^2).^2))','x','y');
d1=fx(Mx,My);
d2=fy(Mx,My);
hold on
subplot(1,1,1)
contour(Mx,My,f,20); colorbar %curvas de nivel función potencial
quiver(Mx,My,d1,d2) %campo de velocidades fluido, ortogonal a curvas de nivel
axis([-5,5,-5,5])
hold off
}}
En nuestro dibujo saldrán las flechas, que representan nuestro campo de velocidades superpuestas con las curvas de nivel de la función potencial. Si ampliamos la imagen podremos observar que estos dos elementos son ortogonales:

[[Archivo: 2campoycurvas.jpg|800px|thumb|center|CAMPO DE VELOCIDADES]]

==º Ortogonalidad <math>\nabla\vec u</math> y <math>\vec n</math> ==
Como el obstáculo es un círculo situado en el plano xy, sabemos con certeza que entonces el vector normal será el correspondiente al eje z, es decir: <math>\vec n= \vec k</math>
Del apartado anterior conocemos que <math>\vec u=\nabla\varphi</math>.

<math>\vec{u}=sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}</math>

Conociendo entonces el vector normal y el vector gradiente, podemos calcular su producto escalar. Para comprobar si dos vectores son ortogonales, su producto escalar debe de ser nulo, es por ello que realizamos el producto escalar del gradiente y del vector normal.
<math>\nabla\varphi\cdot\vec n</math> = <math>(sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}})\cdot\vec ({e}_{z})</math> = 0

Como podemos observar, el producto escalar resulta 0, y queda comprobado así que el vector función potencial es perpendicular a su vector normal.

==º Calculo de <math>\vec u</math> con <math> \frac{1}{\rho}</math> despreciable ==
Al centrarnos en el fluido a una distancia muy lejana, podemos considerar <math> \frac{1}{\rho}</math> despreciable. En este caso sabremos que <math> u</math> cambia. La velocidad de las partículas desde los puntos lejanos del fluido viene definida por el gradiente de la función potencial <math>\varphi=\rho\sin(\theta)</math>

<math>\vec{u}=sin(\theta){\vec{e}_{\rho}}+\frac{1}{\rho}\rho\cos(\theta){\vec{e}_{\theta}}</math> 

En este caso:
<math>|\vec{u}| =\sqrt{cos(\theta)^2+sin(\theta)^2}=1</math>

== º Rotacional y Divergencia ==
:
<math>\nabla\times\vec u=\begin{vmatrix} \vec {e}_{\rho}&\vec {e}_{\theta}&\vec {e}_{z} \\ \frac{\partial}{\partial{\rho}} & \frac{\partial}{\partial{\theta}} & \frac{\partial}{\partial{z}} \\ sin(\theta)[1-\frac{1}{\rho^2}] & cos(\theta)[\rho +\frac{1}{\rho}] & {0} \end{vmatrix}=[-cos(\theta)(1-\frac{1}{\rho^2})]\vec {e}_{z}+[cos(\theta)(1-\frac{1}{\rho^2})]\vec {e}_{z}=\vec {0}</math>
:
<math>\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}({\rho}sin(\theta)(1 -\frac{1}{\rho^2}))+\frac{\partial}{\partial{\theta}}(cos(\theta)(\rho +\frac{1}{\rho})+\frac{\partial}{\partial{z}}(\rho\cdot{0})]=\frac{1}{\rho}[sin(\theta)(1+\frac{1}{\rho^2})-sin(\theta)(1+\frac{1}{\rho^2})]=0 </math>

En este apartado comprobamos analíticamente que la divergencia y rotacional de <math> u</math> son nulas. Que la divergencia sea nula indica que en este caso el volumen del fluido no varía, es decir, se mantiene, ni se expande ni se contrae. A esta condición se le denomina incompresibilidad. Además, que el rotacional sea nulo demuestra que las partículas del fluido no giran.

== ºLíneas de corriente ==
Las líneas de corriente son aquellas que son tangentes a cada punto del campo <math>\vec u</math>. Para ello calcularemos el vector <math>\vec v=\vec k\times\vec u</math>

<math>\vec v=\frac{1}{\rho}\begin{vmatrix} \vec {e}_{\rho}&\rho\vec {e}_{\theta}&\vec {e}_{z} \\ {0} & {0} & {1} \\ sin(\theta)(1-\frac{1}{\rho^2}) & cos(\theta)(1+\frac{1}{\rho^2})\rho & {0} \end{vmatrix}= sin(\theta)(1-\frac{1}{\rho^2})\vec {e}_{\theta} - cos(\theta)(1+\frac{1}{\rho^2})\vec {e}_{\rho} =\vec {grad}\psi</math>:

<math>\psi=\int - cos(\theta)(1+\frac{1}{\rho^2}) \, d\rho \,\!=-cos(\theta)(\rho-\frac{1}{\rho}) + f(\theta) </math>:

<math>\frac{\partial\psi}{\partial\theta}=sin(\theta)(\rho-\frac{1}{\rho})+f'(\theta)=sin(\theta)(\rho-\frac{1}{\rho})\qquad f'(\theta)=0\qquad f(\theta)=cte</math>:

<math>\psi=-cos(\theta)(\rho-\frac{1}{\rho})</math>:

Del apartado anterior, podemos ver que el rotacional es nulo, y por lo tanto es irrotacional. De esta forma, las líneas de corriente serán tangentes al campo. 
Las curvas de nivel de la corriente son las siguientes: 
[[Archivo:curvascorriente.jpg|800px|miniaturadeimagen|centro|LÍNEAS DE CORRIENTE]] 
Si superponemos la imagen con el campo de velocidades antes calculado, observamos que efectivamente son tangentes: 
[[Archivo: curvasycampo.png|800px|miniaturadeimagen|centro|Representación de las líneas de corriente, tangentes al campo <math>\vec u</math>]]

El código Matlab empleado es el siguiente:

{{matlab|codigo=ro=linspace(1,5,300); %intervalo de ro (1,5)
tt=linspace(0,2*pi,300); %intervalo de teta (0,2pi)
[Mro,Mtt]=meshgrid(ro,tt);
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));

syms x y
fi=inline('-x+x./(x.^2+y.^2)','x','y'); %funciÃ³n que define la corriente de u
f2=fi(Mx,My);
fx=inline('-2*(x.*y)./((x.^2+y.^2).^2)','x','y'); %derivadas de funcion potencial que definen campo de velocidades del fluido
fy=inline('1+((x.^2-y.^2)./((x.^2+y.^2).^2))','x','y');
d1=fx(Mx,My);
d2=fy(Mx,My);
hold on
subplot(1,1,1)
contour(Mx,My,f2,100) %curvas de nivel de la corriente
quiver(Mx,My,d1,d2) %campo de velocidades fluido, tangente a lineas de corriente
hold off
}}

== º Puntos frontera==
Al tener definida la superficie del anillo con los radios 1 y 5, buscamos donde están situados los puntos de frontera. Con la función potencial inicial <math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math> y con
<math> u</math> = <math>(sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}) </math>.

Calculamos el modulo del vector en función de <math>\theta</math> para

<math>\rho=1\qquad </math> <math> u</math> = <math> 2cos(\theta){\vec{e}_{\theta}} </math>

<math>|\vec u|=\sqrt{4cos(\theta)^2} = 2cos(\theta) = 2|cos(\theta)|</math> 

Como podemos observar, el módulo de la velocidad es máximo cuando cos<math>\theta</math> alcanza sus valores máximos (en 0 y en <math>\pi</math>, donde alcanza en valor absoluto el valor de 1). 

Así, los puntos de velocidad máxima son <math>\theta=0</math> y <math>\theta={pi}</math>. 

Los puntos de velocidad mínima o puntos de remanso son <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math>. 

== º Presión del fluido ==

En dinámica de fluidos, el principio de Bernoulli, también denominado ecuación de Bernoulli, describe el comportamiento de un fluido moviéndose a lo largo de una línea de corriente. Relacionando las siguientes variables: 

<math>p =</math>presión estática a la que está sometido el fluido debida a las moléculas que lo rodean. 
<math>d =</math>densidad del fluido. 

<math>\vec { u } =</math>velocidad de flujo del fluido. :

<math>\frac { 1 }{ 2 } d { \left| \vec { u } \right| }^{ 2 }+ { p } =cte</math> 
Para que esta expresión pueda aplicarse, el fluido debe ser incompresible, no viscoso y fluir en régimen estacionario (la velocidad en un punto determinado no varía con el tiempo). Tomaremos nuestro fluido como tal. 
Suponiendo que nuestro fluido efectivamente cumple la ecuación de Bernouilli, que posee una densidad constante de d=2 y tomando la cte de la ecuación como cte=10, calcularemos la presión de nuestro fluido. 
Sustituyendo los valores mencionados tenemos que: 
<math> p = 10- |\vec{u}|^2 </math> 
Calculamos el módulo al cuadrado de nuestro campo de velocidades: 
<math> |\vec{u}|^2 = (1-\frac{1}{ρ^2})^2 sen^2(θ) + (1+\frac{1}{ρ^2})^2 cos^2(θ) </math> 
Así la función que define la presión del fluido es: 
<math> p = 10- (1-\frac{1}{ρ^2})^2 sen^2(θ) + (1+\frac{1}{ρ^2})^2 cos^2(θ) </math> 

Para la representación de la presión empleamos el siguiente código:
{{matlab| codigo=
ro=linspace(1,5,50);
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt);
x=Mro.*cos(Mtt);
y=Mro.*sin(Mtt);
P=10-(((1-1./Mro.^2).^2).*(sin(Mtt).^2)+((1+1./Mro.^2).^2).*(cos(Mtt).^2)); %función presión
surf(x,y,P); colorbar;
axis([-5,5,-5,5])
view(2)
Pmax=max(max(P)) %presión máx y mín
Pmin=min(min(P))
}}
Matlab nos devuelve la siguiente imagen y los valores máximos y mínimos de la presión (Presión máx: 9.9959 Presión mínima: 6)
[[Archivo:Presionnn.jpg|800px|miniaturadeimagen|centro|Fig.PRESIÓN DEL FLUIDO]] 
Recordemos ahora la imagen del campo de velocidades (ahora sin ampliar) que habíamos obtenido anteriormente con el código matlab de el apartado 2.2: 

[[Archivo:velocidadd.jpg|800px|miniaturadeimagen|centro|Fig.CAMPO DE VELOCIDADES]] 

Si comparamos nuestra imagen de la presión del fluido (Fig.PRESIÓN DEL FLUIDO) con la del campo de velocidades (Fig.CAMPO DE VELOCIDADES) podemos observar que los puntos que cuentan con módulo de velocidad máxima (<math>\theta=0</math> y <math>\theta={\pi}</math>) son aquellos con presiones mínimas. Esto es un resultado coherente que verifica la ecuación de Bernoulli, a mayor velocidad, menores presiones.

== ºEcuación de Navier-Stokes ==

La ecuación de Navier-Stokes describe el movimiento de los fluidos a partir de su campo de velocidades y de presiones:

<math> d(\vec{u}·\nabla)\vec{u} + \nabla p = \mu \Delta \vec{u} </math> 

Los fluidos incompresibles tienen la propiedad de ocupar siempre el mismo volumen, de esta forma la viscosidad será nula: <math> \mu =0 </math> y la densidad será nula también: <math>d = 2</math>. Si sustituimos en la ecuación general, llegamos a la conclusión de que para un fluido incompresible:

<math> 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 </math> 

Si consideramos los campos <math> (\vec{u},p) </math> hallados anteriormente:

<math> (\vec{u}·\nabla)\vec{u} = u_{1} \frac{\partial u_{1}}{\partial \rho} \vec{e}_{\rho} + u_{1} \frac{\partial u_{1}}{\partial \theta} \vec{e}_{\theta}+ u_{2} \frac{\partial u_{2}}{\partial \rho} \vec{e}_{\rho}+ u_{2} \frac{\partial u_{2}}{\partial \theta} \vec{e}_{\theta} = </math> :

<math> (\vec{u}·\nabla)\vec{u}=(sen(\theta)(1-\frac{1}{\rho^2})[\frac{2}{\rho^3}sin(\theta)\vec{e}_{\rho}+(1-\frac{1}{\rho^2})cos(\theta)\vec{e}_{\theta}])+((1+\frac{1}{\rho^2})cos(\theta)[-\frac{2}{\rho^3}cos(\theta)\vec{e}_{\rho}-(1+\frac{1}{\rho^2})sin(\theta)\vec{e}_{\theta}])=</math>

<math> =(\vec{u}·\nabla)\vec{u}=[\frac{2}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)-\frac{2}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[(1-\frac{1}{\rho^2})^2sin(\theta)cos(\theta)-(1+\frac{1}{\rho^2})^2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

Lo aplicamos a nuestro caso con la densidad determinada:

<math> 2(\vec{u}·\nabla)\vec{u}=[\frac{4}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)- \frac{4}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[(1-\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)-(1+\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

De la misma forma:

<math> \nabla p = \frac{\partial p}{\partial \rho} \vec{e}_{\rho} + \frac{\partial p}{\partial \theta} \vec{e}_{\theta} = </math> :

<math> \nabla p = -[(2(1-\frac{1}{\rho^2})\frac{2}{\rho^3}sin^2(\theta))+(2(1+\frac{1}{\rho^2})(-\frac{2}{\rho^3})cos(\theta))]\vec{e}_{\rho} - [((1-\frac{1}{\rho^2})^2)2sin(\theta)cos(\theta)+((1+\frac{1}{\rho^2})^2)2(-sin(\theta))cos(\theta)]\vec{e}_{\theta}=</math>

<math> = [-\frac{4}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)+ \frac{4}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[-(1-\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)+(1+\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

Habiendo sustituido todos los valores en la ecuación de Navier-Stokes, coparando los resultados vemos directamente que se cumple <math> 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 </math> quedando así demostrado que se cumple la ecuación.

== º Trayectoria de la línea de corriente==
Si fuéramos una partícula del fluido con el que estamos trabajando, seguiríamos la trayectoria de una línea de corriente, es decir empezaríamos con su trayectoria original y a medida que se va acercando al obstáculo recorreríamos el contorno de este, y cuando lo hubiésemos pasado terminado volveríamos a la dirección inicial. En nuestro caso, como el obstáculo tiene un contorno circular el fluido empieza recto, luego se curva para evitar el obstáculo y finalmente vuelve a ir en la dirección inicial.
En relación a la variación de la velocidad y presión podemos comprobar con las gráficas, que en los puntos donde la velocidad es mínima las presiones son máximas. En los puntos que ya hemos calculado anteriormente donde la velocidad es 0, <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math>, sabemos que la presión es máxima.
Con esto podemos deducir que la presión y velocidad de estas partículas son inversamente proporcionales.

== º Paradoja D'Alembert==
:
<math>\vec{u}=sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}</math>

Para <math>\rho=1</math>:

<math>\vec u=2cos(\theta)\vec {e}_{\theta}</math>:

Según el Teorema de Kutta-Joukowski la fuerza que ejerce el fluido sobre el obstáculo, en este caso el círculo de radio 1, es proporcional a la circulación. Por lo tanto:

<math>\int_{0}^{2\pi} \vec u r'(\theta)\,d\theta =\int_{0}^{2\pi} 2cos(\theta)\vec{e}_{\theta}\vec{e}_{\theta}\,d\theta =\int_{0}^{2\pi} 2cos(\theta)\,d\theta=2[sin(0)-sin(2\pi)]=0 </math>:

El fluido no ejerce ninguna fuerza sobre el obstáculo, dando lugar a la paradoja D'Alembert, que concluye que la fuerza resultante sobre el cuerpo es cero, contradiciéndose con la observación. Esto demuestra también que la circulación del campo vectorial a lo largo de la circunferencia se anula.
Para poder ilustrar esta fuerza se necesitan tomar soluciones de modelos de fluidos con una viscosidad positiva.

== º Curvas de nivel de la presión==
Observando las curvas de nivel de la presión podemos entender mejor las conclusiones a las que hemos ido llegando a lo largo del trabajo. En esta gráfica podemos observar los puntos de presión máxima <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math> y los puntos de presión mínima <math>\theta=0</math> y <math>\theta=\pi</math>.
Además, cabe destacar que, en las zonas en las que las curvas de nivel están mas juntas, hay mayor variación de la presión.

[[Archivo:Doce_mnj.jpg|800px|miniaturadeimagen|centro|Líneas de presión]]

Código Matlab:
{{matlab| codigo=

contour(x,y,P,100) % Curvas de nivel de la presión
axis([-5,5,-5,5]) % Región a dibujar

}}
[[Categoría:Teoría de Campos]]
[[Categoría:TC21/22]]

Visualización de campos escalares y vectoriales en fluidos. (Grupo B6)

2021-12-10T11:23:28Z

Julia Cabeza:

{{ TrabajoED | Visualización de campos escalares y vectoriales en fluidos. (Grupo B6) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC21/22|2021-22]] | María Magaldi Jurado, Nuria Martínez Ballester, Julia Cabeza Duque }}

En este trabajo vamos a visualizar campos escalares y vectoriales en fluidos. Más específicamente, el fluido a tratar es uno incompresible que se encuentra alrededor de un obstáculo de forma circular, formando así un anillo comprendido entre las circunferencias de radio 1 y 5, con el centro en el origen de coordenadas. Para facilitar el trabajo usaremos coordenadas polares (cilíndricas <math>(\rho,\theta)</math>). Para dibujar que el fluido va por el exterior de la superficie definida consideraremos que <math>(\rho,\theta) \in [-5,5]*[-5,5]</math>.
Además, tendremos que la velocidad de las partículas del fluido viene definida por el gradiente de la función potencial <math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math>

== º Mallado del fluido ==
El siguiente mallado representa los puntos interiores de la región ocupada por un fluido. Se encuentra en el exterior del círculo unidad. Como ya se ha introducido, nuestro mallado del anillo está comprendido entre los radios 1 y 5 y centro el origen.

[[Archivo:1,mallad0.jpg|800px|thumb|center|Mallado del fluido]]

El CÓDIGO matlab empleado es:
{{matlab|codigo=
ro=linspace(1,5,50); %creación de los intervalos
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt); %creación matrices
Mx=Mro.*cos(Mtt); %parametrización de x
My=Mro.*sin(Mtt); %parametrización de y
Mz=zeros(size(Mx)); %parametrización de z (nula al estar en un plano)
subplot(1,1,1);
mesh(Mx,My,Mz) %mallado
axis([-5,5,-5,5]) %superficie definida en [-5,5]*[-5,5]
view(2)
}}

== º Función potencial <math>\varphi</math> y velocidad de las partículas del fluido ==
=== º Función potencial ===
La velocidad de las partículas del fluido viene dada por el gradiente de la función potencial:
<math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math>
Las curvas de nivel de esta función potencial son las siguientes:
[[Archivo:2curvasdenivel.jpg|800px|thumb|center|FUNCIÓN POTENCIAL]]
El código matlab empleado es:
{{matlab|codigo=
ro=linspace(1,5,50);
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt);
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));
f=(Mro+1./Mro).*sin(Mtt); %función potencial
subplot(1,1,1)
contour(Mx,My,f,50); colorbar %curvas de nivel función potencial
}}

=== º Campo de velocidades (gradiente de función potencial) ===
Si dibujamos nuestro campo de velocidades <math>\vec u=\nabla\varphi</math> en matlab, la representación del mismo con flechas deberá ser ortogonal a las curvas de nivel antes dibujadas. Ya que, de forma geométrica, el gradiente es un vector normal a la curva de nivel de la función de la cual es gradiente en el punto que se esté estudiando.
Hemos pasado nuestra función potencial a cartesianas y trabajado posteriormente con sus derivadas. 
Así, con las igualdades: <math>\rho=\sqrt{x^2+y^2}\qquad\theta=arctg(\frac{y}{x})</math> 
Obtenemos la función potencial: <math>\varphi=y+\frac{y}{x^2+y^2}</math> 

El gradiente se define como:
<math>\vec u=\nabla\varphi=\frac{∂φ}{∂x}\vec i + \frac{∂φ}{∂j}\vec j </math>

Realizamos las derivadas parciales: 
<math>FX=\frac{\partial\varphi(x,y)}{\partial x}=\frac{-2xy}{(x^2+y^2)^2}\qquad FY=\frac{\partial\varphi(x,y)}{\partial y}=\frac{x^2-y^2}{(x^2+y^2)^2}</math> 
Con estas funciones, podremos dibujar en matlab el campo de velocidades con el siguiente código:
{{matlab|codigo=
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));
f=(Mro+1./Mro).*sin(Mtt); %función potencial
syms x y %nuevas variables en cartesianas para definir derivadas (campo de velocidades)
fx=inline('-2*(x.*y)./((x.^2+y.^2).^2)','x','y'); %derivadas que definen nuestro campos de velocidades
fy=inline('1+((x.^2-y.^2)./((x.^2+y.^2).^2))','x','y');
d1=fx(Mx,My);
d2=fy(Mx,My);
hold on
subplot(1,1,1)
contour(Mx,My,f,20); colorbar %curvas de nivel función potencial
quiver(Mx,My,d1,d2) %campo de velocidades fluido, ortogonal a curvas de nivel
axis([-5,5,-5,5])
hold off
}}
En nuestro dibujo saldrán las flechas, que representan nuestro campo de velocidades superpuestas con las curvas de nivel de la función potencial. Si ampliamos la imagen podremos observar que estos dos elementos son ortogonales:

[[Archivo: 2campoycurvas.jpg|800px|thumb|center|CAMPO DE VELOCIDADES]]

==º Ortogonalidad <math>\nabla\vec u</math> y <math>\vec n</math> ==
Como el obstáculo es un círculo situado en el plano xy, sabemos con certeza que entonces el vector normal será el correspondiente al eje z, es decir: <math>\vec n= \vec k</math>
Del apartado anterior conocemos que <math>\vec u=\nabla\varphi</math>.

<math>\vec{u}=sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}</math>

Conociendo entonces el vector normal y el vector gradiente, podemos calcular su producto escalar. Para comprobar si dos vectores son ortogonales, su producto escalar debe de ser nulo, es por ello que realizamos el producto escalar del gradiente y del vector normal.
<math>\nabla\varphi\cdot\vec n</math> = <math>(sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}})\cdot\vec ({e}_{z})</math> = 0

Como podemos observar, el producto escalar resulta 0, y queda comprobado así que el vector función potencial es perpendicular a su vector normal.

==º Calculo de <math>\vec u</math> con <math> \frac{1}{\rho}</math> despreciable ==
Al centrarnos en el fluido a una distancia muy lejana, podemos considerar <math> \frac{1}{\rho}</math> despreciable. En este caso sabremos que <math> u</math> cambia. La velocidad de las partículas desde los puntos lejanos del fluido viene definida por el gradiente de la función potencial <math>\varphi=\rho\sin(\theta)</math>

<math>\vec{u}=sin(\theta){\vec{e}_{\rho}}+\frac{1}{\rho}\rho\cos(\theta){\vec{e}_{\theta}}</math> 

En este caso:
<math>|\vec{u}| =\sqrt{cos(\theta)^2+sin(\theta)^2}=1</math>

== º Rotacional y Divergencia ==
:
<math>\nabla\times\vec u=\begin{vmatrix} \vec {e}_{\rho}&\vec {e}_{\theta}&\vec {e}_{z} \\ \frac{\partial}{\partial{\rho}} & \frac{\partial}{\partial{\theta}} & \frac{\partial}{\partial{z}} \\ sin(\theta)[1-\frac{1}{\rho^2}] & cos(\theta)[\rho +\frac{1}{\rho}] & {0} \end{vmatrix}=[-cos(\theta)(1-\frac{1}{\rho^2})]\vec {e}_{z}+[cos(\theta)(1-\frac{1}{\rho^2})]\vec {e}_{z}=\vec {0}</math>
:
<math>\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}({\rho}sin(\theta)(1 -\frac{1}{\rho^2}))+\frac{\partial}{\partial{\theta}}(cos(\theta)(\rho +\frac{1}{\rho})+\frac{\partial}{\partial{z}}(\rho\cdot{0})]=\frac{1}{\rho}[sin(\theta)(1+\frac{1}{\rho^2})-sin(\theta)(1+\frac{1}{\rho^2})]=0 </math>

En este apartado comprobamos analíticamente que la divergencia y rotacional de <math> u</math> son nulas. Que la divergencia sea nula indica que en este caso el volumen del fluido no varía, es decir, se mantiene, ni se expande ni se contrae. A esta condición se le denomina incompresibilidad. Además, que el rotacional sea nulo demuestra que las partículas del fluido no giran.

== ºLíneas de corriente ==
Las líneas de corriente son aquellas que son tangentes a cada punto del campo <math>\vec u</math>. Para ello calcularemos el vector <math>\vec v=\vec k\times\vec u</math>

<math>\vec v=\frac{1}{\rho}\begin{vmatrix} \vec {e}_{\rho}&\rho\vec {e}_{\theta}&\vec {e}_{z} \\ {0} & {0} & {1} \\ sin(\theta)(1-\frac{1}{\rho^2}) & cos(\theta)(1+\frac{1}{\rho^2})\rho & {0} \end{vmatrix}= sin(\theta)(1-\frac{1}{\rho^2})\vec {e}_{\theta} - cos(\theta)(1+\frac{1}{\rho^2})\vec {e}_{\rho} =\vec {grad}\psi</math>:

<math>\psi=\int - cos(\theta)(1+\frac{1}{\rho^2}) \, d\rho \,\!=-cos(\theta)(\rho-\frac{1}{\rho}) + f(\theta) </math>:

<math>\frac{\partial\psi}{\partial\theta}=sin(\theta)(\rho-\frac{1}{\rho})+f'(\theta)=sin(\theta)(\rho-\frac{1}{\rho})\qquad f'(\theta)=0\qquad f(\theta)=cte</math>:

<math>\psi=-cos(\theta)(\rho-\frac{1}{\rho})</math>:

Del apartado anterior, podemos ver que el rotacional es nulo, y por lo tanto es irrotacional. De esta forma, las líneas de corriente serán tangentes al campo. 
Las curvas de nivel de la corriente son las siguientes: 
[[Archivo:curvascorriente.jpg|800px|miniaturadeimagen|centro|LÍNEAS DE CORRIENTE]] 
Si superponemos la imagen con el campo de velocidades antes calculado, observamos que efectivamente son tangentes: 
[[Archivo: curvasycampo.png|800px|miniaturadeimagen|centro|Representación de las líneas de corriente, tangentes al campo <math>\vec u</math>]]

El código Matlab empleado es el siguiente:

{{matlab|codigo=ro=linspace(1,5,300); %intervalo de ro (1,5)
tt=linspace(0,2*pi,300); %intervalo de teta (0,2pi)
[Mro,Mtt]=meshgrid(ro,tt);
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));

syms x y
fi=inline('-x+x./(x.^2+y.^2)','x','y'); %funciÃ³n que define la corriente de u
f2=fi(Mx,My);
fx=inline('-2*(x.*y)./((x.^2+y.^2).^2)','x','y'); %derivadas de funcion potencial que definen campo de velocidades del fluido
fy=inline('1+((x.^2-y.^2)./((x.^2+y.^2).^2))','x','y');
d1=fx(Mx,My);
d2=fy(Mx,My);
hold on
subplot(1,1,1)
contour(Mx,My,f2,100) %curvas de nivel de la corriente
quiver(Mx,My,d1,d2) %campo de velocidades fluido, tangente a lineas de corriente
hold off
}}

== º Puntos frontera==
Al tener definida la superficie del anillo con los radios 1 y 5, buscamos donde están situados los puntos de frontera. Con la función potencial inicial <math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math> y con
<math> u</math> = <math>(sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}) </math>.

Calculamos el modulo del vector en función de <math>\theta</math> para

<math>\rho=1\qquad </math> <math> u</math> = <math> 2cos(\theta){\vec{e}_{\theta}} </math>

<math>|\vec u|=\sqrt{4cos(\theta)^2} = 2cos(\theta) = 2|cos(\theta)|</math> 

Como podemos observar, el módulo de la velocidad es máximo cuando cos<math>\theta</math> alcanza sus valores máximos (en 0 y en <math>\pi</math>, donde alcanza en valor absoluto el valor de 1). 

Así, los puntos de velocidad máxima son <math>\theta=0</math> y <math>\theta={pi}</math>. 

Los puntos de velocidad mínima o puntos de remanso son <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math>. 

== º Presión del fluido ==

En dinámica de fluidos, el principio de Bernoulli, también denominado ecuación de Bernoulli, describe el comportamiento de un fluido moviéndose a lo largo de una línea de corriente. Relacionando las siguientes variables: 

<math>p =</math>presión estática a la que está sometido el fluido debida a las moléculas que lo rodean. 
<math>d =</math>densidad del fluido. 

<math>\vec { u } =</math>velocidad de flujo del fluido. :

<math>\frac { 1 }{ 2 } d { \left| \vec { u } \right| }^{ 2 }+ { p } =cte</math> 
Para que esta expresión pueda aplicarse, el fluido debe ser incompresible, no viscoso y fluir en régimen estacionario (la velocidad en un punto determinado no varía con el tiempo). Tomaremos nuestro fluido como tal. 
Suponiendo que nuestro fluido efectivamente cumple la ecuación de Bernouilli, que posee una densidad constante de d=2 y tomando la cte de la ecuación como cte=10, calcularemos la presión de nuestro fluido. 
Sustituyendo los valores mencionados tenemos que: 
<math> p = 10- |\vec{u}|^2 </math> 
Calculamos el módulo al cuadrado de nuestro campo de velocidades: 
<math> |\vec{u}|^2 = (1-\frac{1}{ρ^2})^2 sen^2(θ) + (1+\frac{1}{ρ^2})^2 cos^2(θ) </math> 
Así la función que define la presión del fluido es: 
<math> p = 10- (1-\frac{1}{ρ^2})^2 sen^2(θ) + (1+\frac{1}{ρ^2})^2 cos^2(θ) </math> 

Para la representación de la presión empleamos el siguiente código:
{{matlab| codigo=
ro=linspace(1,5,50);
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt);
x=Mro.*cos(Mtt);
y=Mro.*sin(Mtt);
P=10-(((1-1./Mro.^2).^2).*(sin(Mtt).^2)+((1+1./Mro.^2).^2).*(cos(Mtt).^2)); %función presión
surf(x,y,P); colorbar;
axis([-5,5,-5,5])
view(2)
Pmax=max(max(P)) %presión máx y mín
Pmin=min(min(P))
}}
Matlab nos devuelve la siguiente imagen y los valores máximos y mínimos de la presión (Presión máx: 9.9959 Presión mínima: 6)
[[Archivo:Presionnn.jpg|800px|miniaturadeimagen|centro|Fig.PRESIÓN DEL FLUIDO]] 
Recordemos ahora la imagen del campo de velocidades (ahora sin ampliar) que habíamos obtenido anteriormente con el código matlab de el apartado 2.2: 

[[Archivo:velocidadd.jpg|800px|miniaturadeimagen|centro|Fig.CAMPO DE VELOCIDADES]] 

Si comparamos nuestra imagen de la presión del fluido (Fig.PRESIÓN DEL FLUIDO) con la del campo de velocidades (Fig.CAMPO DE VELOCIDADES) podemos observar que los puntos que cuentan con módulo de velocidad máxima (<math>\theta=0</math> y <math>\theta={\pi}</math>) son aquellos con presiones mínimas. Esto es un resultado coherente que verifica la ecuación de Bernoulli, a mayor velocidad, menores presiones.

== ºEcuación de Navier-Stokes ==

La ecuación de Navier-Stokes describe el movimiento de los fluidos a partir de su campo de velocidades y de presiones:

<math> d(\vec{u}·\nabla)\vec{u} + \nabla p = \mu \Delta \vec{u} </math> 

Los fluidos incompresibles tienen la propiedad de ocupar siempre el mismo volumen, de esta forma la viscosidad será nula: <math> \mu =0 </math> y la densidad será nula también: <math>d = 2</math>. Si sustituimos en la ecuación general, llegamos a la conclusión de que para un fluido incompresible:

<math> 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 </math> 

Si consideramos los campos <math> (\vec{u},p) </math> hallados anteriormente:

<math> (\vec{u}·\nabla)\vec{u} = u_{1} \frac{\partial u_{1}}{\partial \rho} \vec{e}_{\rho} + u_{1} \frac{\partial u_{1}}{\partial \theta} \vec{e}_{\theta}+ u_{2} \frac{\partial u_{2}}{\partial \rho} \vec{e}_{\rho}+ u_{2} \frac{\partial u_{2}}{\partial \theta} \vec{e}_{\theta} = </math> :

<math> (\vec{u}·\nabla)\vec{u}=(sen(\theta)(1-\frac{1}{\rho^2})[\frac{2}{\rho^3}sin(\theta)\vec{e}_{\rho}+(1-\frac{1}{\rho^2})cos(\theta)\vec{e}_{\theta}])+((1+\frac{1}{\rho^2})cos(\theta)[-\frac{2}{\rho^3}cos(\theta)\vec{e}_{\rho}-(1+\frac{1}{\rho^2})sin(\theta)\vec{e}_{\theta}])=</math>

<math> =(\vec{u}·\nabla)\vec{u}=[\frac{2}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)-\frac{2}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[(1-\frac{1}{\rho^2})^2sin(\theta)cos(\theta)-(1+\frac{1}{\rho^2})^2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

Lo aplicamos a nuestro caso con la densidad determinada:

<math> 2(\vec{u}·\nabla)\vec{u}=[\frac{4}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)- \frac{4}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[(1-\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)-(1+\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

De la misma forma:

<math> \nabla p = \frac{\partial p}{\partial \rho} \vec{e}_{\rho} + \frac{\partial p}{\partial \theta} \vec{e}_{\theta} = </math> :

<math> \nabla p = -[(2(1-\frac{1}{\rho^2})\frac{2}{\rho^3}sin^2(\theta))+(2(1+\frac{1}{\rho^2})(-\frac{2}{\rho^3})cos(\theta))]\vec{e}_{\rho} - [((1-\frac{1}{\rho^2})^2)2sin(\theta)cos(\theta)+((1+\frac{1}{\rho^2})^2)2(-sin(\theta))cos(\theta)]\vec{e}_{\theta}=</math>

<math> = [-\frac{4}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)+ \frac{4}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[-(1-\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)+(1+\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

Habiendo sustituido todos los valores en la ecuación de Navier-Stokes, coparando los resultados vemos directamente que se cumple <math> 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 </math> quedando así demostrado que se cumple la ecuación.

== º Trayectoria de la línea de corriente==
Si fuéramos una partícula del fluido con el que estamos trabajando, seguiríamos la trayectoria de una línea de corriente, es decir empezaríamos con su trayectoria original y a medida que se va acercando al obstáculo recorreríamos el contorno de este, y cuando lo hubiésemos pasado terminado volveríamos a la dirección inicial. En nuestro caso, como el obstáculo tiene un contorno circular el fluido empieza recto, luego se curva para evitar el obstáculo y finalmente vuelve a ir en la dirección inicial.
En relación a la variación de la velocidad y presión podemos comprobar con las gráficas, que en los puntos donde la velocidad es mínima las presiones son máximas. En los puntos que ya hemos calculado anteriormente donde la velocidad es 0, <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math>, sabemos que la presión es máxima.
Con esto podemos deducir que la presión y velocidad de estas partículas son inversamente proporcionales.

== º Paradoja D'Alembert==
:
<math>\vec{u}=sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}</math>

Para <math>\rho=1</math>:

<math>\vec u=2cos(\theta)\vec {e}_{\theta}</math>:

Según el Teorema de Kutta-Joukowski la fuerza que ejerce el fluido sobre el obstáculo, en este caso el círculo de radio 1, es proporcional a la circulación. Por lo tanto:

<math>\int_{0}^{2\pi} \vec u r'(\theta)\,d\theta =\int_{0}^{2\pi} 2cos(\theta)\vec{e}_{\theta}\vec{e}_{\theta}\,d\theta =\int_{0}^{2\pi} 2cos(\theta)\,d\theta=2[sin(0)-sin(2\pi)]=0 </math>:

El fluido no ejerce ninguna fuerza sobre el obstáculo, dando lugar a la paradoja D'Alembert, que concluye que la fuerza resultante sobre el cuerpo es cero, contradiciéndose con la observación. Esto demuestra también que la circulación del campo vectorial a lo largo de la circunferencia se anula.
Para poder ilustrar esta fuerza se necesitan tomar soluciones de modelos de fluidos con una viscosidad positiva.

== º Curvas de nivel de la presión==
Observando las curvas de nivel de la presión podemos entender mejor las conclusiones a las que hemos ido llegando a lo largo del trabajo. En esta gráfica podemos observar los puntos de presión máxima <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math> y los puntos de presión mínima <math>\theta=0</math> y <math>\theta=\pi</math>.
Además, cabe destacar que, en las zonas en las que las curvas de nivel están mas juntas, hay mayor variación de la presión.

[[Archivo:Doce_mnj.jpg|800px|miniaturadeimagen|centro|Líneas de presión]]

Código Matlab:
{{matlab| codigo=

contour(x,y,P,100) % Curvas de nivel de la presión
axis([-5,5,-5,5]) % Región a dibujar

}}
[[Categoría:Teoría de Campos]]
[[Categoría:TC21/22]]

Archivo:Curvascorriente.jpg

2021-12-10T11:19:03Z

Julia Cabeza:

Visualización de campos escalares y vectoriales en fluidos. (Grupo B6)

2021-12-10T10:51:08Z

Julia Cabeza:

{{ TrabajoED | Visualización de campos escalares y vectoriales en fluidos. (Grupo B6) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC21/22|2021-22]] | María Magaldi Jurado, Nuria Martínez Ballester, Julia Cabeza Duque }}

En este trabajo vamos a visualizar campos escalares y vectoriales en fluidos. Más específicamente, el fluido a tratar es uno incompresible que se encuentra alrededor de un obstáculo de forma circular, formando así un anillo comprendido entre las circunferencias de radio 1 y 5, con el centro en el origen de coordenadas. Para facilitar el trabajo usaremos coordenadas polares (cilíndricas <math>(\rho,\theta)</math>). Para dibujar que el fluido va por el exterior de la superficie definida consideraremos que <math>(\rho,\theta) \in [-5,5]*[-5,5]</math>.
Además, tendremos que la velocidad de las partículas del fluido viene definida por el gradiente de la función potencial <math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math>

== º Mallado del fluido ==
El siguiente mallado representa los puntos interiores de la región ocupada por un fluido. Se encuentra en el exterior del círculo unidad. Como ya se ha introducido, nuestro mallado del anillo está comprendido entre los radios 1 y 5 y centro el origen.

[[Archivo:1,mallad0.jpg|800px|thumb|center|Mallado del fluido]]

El CÓDIGO matlab empleado es:
{{matlab|codigo=
ro=linspace(1,5,50); %creación de los intervalos
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt); %creación matrices
Mx=Mro.*cos(Mtt); %parametrización de x
My=Mro.*sin(Mtt); %parametrización de y
Mz=zeros(size(Mx)); %parametrización de z (nula al estar en un plano)
subplot(1,1,1);
mesh(Mx,My,Mz) %mallado
axis([-5,5,-5,5]) %superficie definida en [-5,5]*[-5,5]
view(2)
}}

== º Función potencial <math>\varphi</math> y velocidad de las partículas del fluido ==
=== º Función potencial ===
La velocidad de las partículas del fluido viene dada por el gradiente de la función potencial:
<math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math>
Las curvas de nivel de esta función potencial son las siguientes:
[[Archivo:2curvasdenivel.jpg|800px|thumb|center|FUNCIÓN POTENCIAL]]
El código matlab empleado es:
{{matlab|codigo=
ro=linspace(1,5,50);
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt);
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));
f=(Mro+1./Mro).*sin(Mtt); %función potencial
subplot(1,1,1)
contour(Mx,My,f,50); colorbar %curvas de nivel función potencial
}}

=== º Campo de velocidades (gradiente de función potencial) ===
Si dibujamos nuestro campo de velocidades <math>\vec u=\nabla\varphi</math> en matlab, la representación del mismo con flechas deberá ser ortogonal a las curvas de nivel antes dibujadas. Ya que, de forma geométrica, el gradiente es un vector normal a la curva de nivel de la función de la cual es gradiente en el punto que se esté estudiando.
Hemos pasado nuestra función potencial a cartesianas y trabajado posteriormente con sus derivadas. 
Así, con las igualdades: <math>\rho=\sqrt{x^2+y^2}\qquad\theta=arctg(\frac{y}{x})</math> 
Obtenemos la función potencial: <math>\varphi=y+\frac{y}{x^2+y^2}</math> 

El gradiente se define como:
<math>\vec u=\nabla\varphi=\frac{∂φ}{∂x}\vec i + \frac{∂φ}{∂j}\vec j </math>

Realizamos las derivadas parciales: 
<math>FX=\frac{\partial\varphi(x,y)}{\partial x}=\frac{-2xy}{(x^2+y^2)^2}\qquad FY=\frac{\partial\varphi(x,y)}{\partial y}=\frac{x^2-y^2}{(x^2+y^2)^2}</math> 
Con estas funciones, podremos dibujar en matlab el campo de velocidades con el siguiente código:
{{matlab|codigo=
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));
f=(Mro+1./Mro).*sin(Mtt); %función potencial
syms x y %nuevas variables en cartesianas para definir derivadas (campo de velocidades)
fx=inline('-2*(x.*y)./((x.^2+y.^2).^2)','x','y'); %derivadas que definen nuestro campos de velocidades
fy=inline('1+((x.^2-y.^2)./((x.^2+y.^2).^2))','x','y');
d1=fx(Mx,My);
d2=fy(Mx,My);
hold on
subplot(1,1,1)
contour(Mx,My,f,20); colorbar %curvas de nivel función potencial
quiver(Mx,My,d1,d2) %campo de velocidades fluido, ortogonal a curvas de nivel
axis([-5,5,-5,5])
hold off
}}
En nuestro dibujo saldrán las flechas, que representan nuestro campo de velocidades superpuestas con las curvas de nivel de la función potencial. Si ampliamos la imagen podremos observar que estos dos elementos son ortogonales:

[[Archivo: 2campoycurvas.jpg|800px|thumb|center|CAMPO DE VELOCIDADES]]

==º Ortogonalidad <math>\nabla\vec u</math> y <math>\vec n</math> ==
Como el obstáculo es un círculo situado en el plano xy, sabemos con certeza que entonces el vector normal será el correspondiente al eje z, es decir: <math>\vec n= \vec k</math>
Del apartado anterior conocemos que <math>\vec u=\nabla\varphi</math>.

<math>\vec{u}=sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}</math>

Conociendo entonces el vector normal y el vector gradiente, podemos calcular su producto escalar. Para comprobar si dos vectores son ortogonales, su producto escalar debe de ser nulo, es por ello que realizamos el producto escalar del gradiente y del vector normal.
<math>\nabla\varphi\cdot\vec n</math> = <math>(sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}})\cdot\vec ({e}_{z})</math> = 0

Como podemos observar, el producto escalar resulta 0, y queda comprobado así que el vector función potencial es perpendicular a su vector normal.

==º Calculo de <math>\vec u</math> con <math> \frac{1}{\rho}</math> despreciable ==
Al centrarnos en el fluido a una distancia muy lejana, podemos considerar <math> \frac{1}{\rho}</math> despreciable. En este caso sabremos que <math> u</math> cambia. La velocidad de las partículas desde los puntos lejanos del fluido viene definida por el gradiente de la función potencial <math>\varphi=\rho\sin(\theta)</math>

<math>\vec{u}=sin(\theta){\vec{e}_{\rho}}+\frac{1}{\rho}\rho\cos(\theta){\vec{e}_{\theta}}</math> 

En este caso:
<math>|\vec{u}| =\sqrt{cos(\theta)^2+sin(\theta)^2}=1</math>

== º Rotacional y Divergencia ==
:
<math>\nabla\times\vec u=\begin{vmatrix} \vec {e}_{\rho}&\vec {e}_{\theta}&\vec {e}_{z} \\ \frac{\partial}{\partial{\rho}} & \frac{\partial}{\partial{\theta}} & \frac{\partial}{\partial{z}} \\ sin(\theta)[1-\frac{1}{\rho^2}] & cos(\theta)[\rho +\frac{1}{\rho}] & {0} \end{vmatrix}=[-cos(\theta)(1-\frac{1}{\rho^2})]\vec {e}_{z}+[cos(\theta)(1-\frac{1}{\rho^2})]\vec {e}_{z}=\vec {0}</math>
:
<math>\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}({\rho}sin(\theta)(1 -\frac{1}{\rho^2}))+\frac{\partial}{\partial{\theta}}(cos(\theta)(\rho +\frac{1}{\rho})+\frac{\partial}{\partial{z}}(\rho\cdot{0})]=\frac{1}{\rho}[sin(\theta)(1+\frac{1}{\rho^2})-sin(\theta)(1+\frac{1}{\rho^2})]=0 </math>

En este apartado comprobamos analíticamente que la divergencia y rotacional de <math> u</math> son nulas. Que la divergencia sea nula indica que en este caso el volumen del fluido no varía, es decir, se mantiene, ni se expande ni se contrae. A esta condición se le denomina incompresibilidad. Además, que el rotacional sea nulo demuestra que las partículas del fluido no giran.

== ºLíneas de corriente ==
Las líneas de corriente son aquellas que son tangentes a cada punto del campo <math>\vec u</math>. Para ello calcularemos el vector <math>\vec v=\vec k\times\vec u</math>

<math>\vec v=\frac{1}{\rho}\begin{vmatrix} \vec {e}_{\rho}&\rho\vec {e}_{\theta}&\vec {e}_{z} \\ {0} & {0} & {1} \\ sin(\theta)(1-\frac{1}{\rho^2}) & cos(\theta)(1+\frac{1}{\rho^2})\rho & {0} \end{vmatrix}= sin(\theta)(1-\frac{1}{\rho^2})\vec {e}_{\theta} - cos(\theta)(1+\frac{1}{\rho^2})\vec {e}_{\rho} =\vec {grad}\psi</math>:

<math>\psi=\int - cos(\theta)(1+\frac{1}{\rho^2}) \, d\rho \,\!=-cos(\theta)(\rho-\frac{1}{\rho}) + f(\theta) </math>:

<math>\frac{\partial\psi}{\partial\theta}=sin(\theta)(\rho-\frac{1}{\rho})+f'(\theta)=sin(\theta)(\rho-\frac{1}{\rho})\qquad f'(\theta)=0\qquad f(\theta)=cte</math>:

<math>\psi=-cos(\theta)(\rho-\frac{1}{\rho})</math>:

Del apartado anterior, podemos ver que el rotacional es nulo, y por lo tanto es irrotacional. De esta forma, las líneas de corriente serán tangentes al campo.

[[Archivo: curvasycampo.png|800px|miniaturadeimagen|centro|Representación de las líneas de corriente, tangentes al campo <math>\vec u</math>]]

El código Matlab empleado es el siguiente:

{{matlab|codigo=ro=linspace(1,5,300); %intervalo de ro (1,5)
tt=linspace(0,2*pi,300); %intervalo de teta (0,2pi)
[Mro,Mtt]=meshgrid(ro,tt);
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));

syms x y
fi=inline('-x+x./(x.^2+y.^2)','x','y'); %funciÃ³n que define la corriente de u
f2=fi(Mx,My);
fx=inline('-2*(x.*y)./((x.^2+y.^2).^2)','x','y'); %derivadas de funcion potencial que definen campo de velocidades del fluido
fy=inline('1+((x.^2-y.^2)./((x.^2+y.^2).^2))','x','y');
d1=fx(Mx,My);
d2=fy(Mx,My);
hold on
subplot(1,1,1)
contour(Mx,My,f2,100) %curvas de nivel de la corriente
quiver(Mx,My,d1,d2) %campo de velocidades fluido, tangente a lineas de corriente
hold off
}}

== º Puntos frontera==
Al tener definida la superficie del anillo con los radios 1 y 5, buscamos donde están situados los puntos de frontera. Con la función potencial inicial <math>\varphi=(\rho+\frac{1}{\rho})sin(\theta)</math> y con
<math> u</math> = <math>(sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}) </math>.

Calculamos el modulo del vector en función de <math>\theta</math> para

<math>\rho=1\qquad </math> <math> u</math> = <math> 2cos(\theta){\vec{e}_{\theta}} </math>

<math>|\vec u|=\sqrt{4cos(\theta)^2} = 2cos(\theta) = 2|cos(\theta)|</math> 

Como podemos observar, el módulo de la velocidad es máximo cuando cos<math>\theta</math> alcanza sus valores máximos (en 0 y en <math>\pi</math>, donde alcanza en valor absoluto el valor de 1). 

Así, los puntos de velocidad máxima son <math>\theta=0</math> y <math>\theta={pi}</math>. 

Los puntos de velocidad mínima o puntos de remanso son <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math>. 

== º Presión del fluido ==

En dinámica de fluidos, el principio de Bernoulli, también denominado ecuación de Bernoulli, describe el comportamiento de un fluido moviéndose a lo largo de una línea de corriente. Relacionando las siguientes variables: 

<math>p =</math>presión estática a la que está sometido el fluido debida a las moléculas que lo rodean. 
<math>d =</math>densidad del fluido. 

<math>\vec { u } =</math>velocidad de flujo del fluido. :

<math>\frac { 1 }{ 2 } d { \left| \vec { u } \right| }^{ 2 }+ { p } =cte</math> 
Para que esta expresión pueda aplicarse, el fluido debe ser incompresible, no viscoso y fluir en régimen estacionario (la velocidad en un punto determinado no varía con el tiempo). Tomaremos nuestro fluido como tal. 
Suponiendo que nuestro fluido efectivamente cumple la ecuación de Bernouilli, que posee una densidad constante de d=2 y tomando la cte de la ecuación como cte=10, calcularemos la presión de nuestro fluido. 
Sustituyendo los valores mencionados tenemos que: 
<math> p = 10- |\vec{u}|^2 </math> 
Calculamos el módulo al cuadrado de nuestro campo de velocidades: 
<math> |\vec{u}|^2 = (1-\frac{1}{ρ^2})^2 sen^2(θ) + (1+\frac{1}{ρ^2})^2 cos^2(θ) </math> 
Así la función que define la presión del fluido es: 
<math> p = 10- (1-\frac{1}{ρ^2})^2 sen^2(θ) + (1+\frac{1}{ρ^2})^2 cos^2(θ) </math> 

Para la representación de la presión empleamos el siguiente código:
{{matlab| codigo=
ro=linspace(1,5,50);
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt);
x=Mro.*cos(Mtt);
y=Mro.*sin(Mtt);
P=10-(((1-1./Mro.^2).^2).*(sin(Mtt).^2)+((1+1./Mro.^2).^2).*(cos(Mtt).^2)); %función presión
surf(x,y,P); colorbar;
axis([-5,5,-5,5])
view(2)
Pmax=max(max(P)) %presión máx y mín
Pmin=min(min(P))
}}
Matlab nos devuelve la siguiente imagen y los valores máximos y mínimos de la presión (Presión máx: 9.9959 Presión mínima: 6)
[[Archivo:Presionnn.jpg|800px|miniaturadeimagen|centro|Fig.PRESIÓN DEL FLUIDO]] 
Recordemos ahora la imagen del campo de velocidades (ahora sin ampliar) que habíamos obtenido anteriormente con el código matlab de el apartado 2.2: 

[[Archivo:velocidadd.jpg|800px|miniaturadeimagen|centro|Fig.CAMPO DE VELOCIDADES]] 

Si comparamos nuestra imagen de la presión del fluido (Fig.PRESIÓN DEL FLUIDO) con la del campo de velocidades (Fig.CAMPO DE VELOCIDADES) podemos observar que los puntos que cuentan con módulo de velocidad máxima (<math>\theta=0</math> y <math>\theta={\pi}</math>) son aquellos con presiones mínimas. Esto es un resultado coherente que verifica la ecuación de Bernoulli, a mayor velocidad, menores presiones.

== ºEcuación de Navier-Stokes ==

La ecuación de Navier-Stokes describe el movimiento de los fluidos a partir de su campo de velocidades y de presiones:

<math> d(\vec{u}·\nabla)\vec{u} + \nabla p = \mu \Delta \vec{u} </math> 

Los fluidos incompresibles tienen la propiedad de ocupar siempre el mismo volumen, de esta forma la viscosidad será nula: <math> \mu =0 </math> y la densidad será nula también: <math>d = 2</math>. Si sustituimos en la ecuación general, llegamos a la conclusión de que para un fluido incompresible:

<math> 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 </math> 

Si consideramos los campos <math> (\vec{u},p) </math> hallados anteriormente:

<math> (\vec{u}·\nabla)\vec{u} = u_{1} \frac{\partial u_{1}}{\partial \rho} \vec{e}_{\rho} + u_{1} \frac{\partial u_{1}}{\partial \theta} \vec{e}_{\theta}+ u_{2} \frac{\partial u_{2}}{\partial \rho} \vec{e}_{\rho}+ u_{2} \frac{\partial u_{2}}{\partial \theta} \vec{e}_{\theta} = </math> :

<math> (\vec{u}·\nabla)\vec{u}=(sen(\theta)(1-\frac{1}{\rho^2})[\frac{2}{\rho^3}sin(\theta)\vec{e}_{\rho}+(1-\frac{1}{\rho^2})cos(\theta)\vec{e}_{\theta}])+((1+\frac{1}{\rho^2})cos(\theta)[-\frac{2}{\rho^3}cos(\theta)\vec{e}_{\rho}-(1+\frac{1}{\rho^2})sin(\theta)\vec{e}_{\theta}])=</math>

<math> =(\vec{u}·\nabla)\vec{u}=[\frac{2}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)-\frac{2}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[(1-\frac{1}{\rho^2})^2sin(\theta)cos(\theta)-(1+\frac{1}{\rho^2})^2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

Lo aplicamos a nuestro caso con la densidad determinada:

<math> 2(\vec{u}·\nabla)\vec{u}=[\frac{4}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)- \frac{4}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[(1-\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)-(1+\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

De la misma forma:

<math> \nabla p = \frac{\partial p}{\partial \rho} \vec{e}_{\rho} + \frac{\partial p}{\partial \theta} \vec{e}_{\theta} = </math> :

<math> \nabla p = -[(2(1-\frac{1}{\rho^2})\frac{2}{\rho^3}sin^2(\theta))+(2(1+\frac{1}{\rho^2})(-\frac{2}{\rho^3})cos(\theta))]\vec{e}_{\rho} - [((1-\frac{1}{\rho^2})^2)2sin(\theta)cos(\theta)+((1+\frac{1}{\rho^2})^2)2(-sin(\theta))cos(\theta)]\vec{e}_{\theta}=</math>

<math> = [-\frac{4}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)+ \frac{4}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[-(1-\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)+(1+\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)]\vec{e}_{\theta} </math> 

Habiendo sustituido todos los valores en la ecuación de Navier-Stokes, coparando los resultados vemos directamente que se cumple <math> 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 </math> quedando así demostrado que se cumple la ecuación.

== º Trayectoria de la línea de corriente==
Si fuéramos una partícula del fluido con el que estamos trabajando, seguiríamos la trayectoria de una línea de corriente, es decir empezaríamos con su trayectoria original y a medida que se va acercando al obstáculo recorreríamos el contorno de este, y cuando lo hubiésemos pasado terminado volveríamos a la dirección inicial. En nuestro caso, como el obstáculo tiene un contorno circular el fluido empieza recto, luego se curva para evitar el obstáculo y finalmente vuelve a ir en la dirección inicial.
En relación a la variación de la velocidad y presión podemos comprobar con las gráficas, que en los puntos donde la velocidad es mínima las presiones son máximas. En los puntos que ya hemos calculado anteriormente donde la velocidad es 0, <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math>, sabemos que la presión es máxima.
Con esto podemos deducir que la presión y velocidad de estas partículas son inversamente proporcionales.

== º Paradoja D'Alembert==
:
<math>\vec{u}=sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}</math>

Para <math>\rho=1</math>:

<math>\vec u=2cos(\theta)\vec {e}_{\theta}</math>:

Según el Teorema de Kutta-Joukowski la fuerza que ejerce el fluido sobre el obstáculo, en este caso el círculo de radio 1, es proporcional a la circulación. Por lo tanto:

<math>\int_{0}^{2\pi} \vec u r'(\theta)\,d\theta =\int_{0}^{2\pi} 2cos(\theta)\vec{e}_{\theta}\vec{e}_{\theta}\,d\theta =\int_{0}^{2\pi} 2cos(\theta)\,d\theta=2[sin(0)-sin(2\pi)]=0 </math>:

El fluido no ejerce ninguna fuerza sobre el obstáculo, dando lugar a la paradoja D'Alembert, que concluye que la fuerza resultante sobre el cuerpo es cero, contradiciéndose con la observación. Esto demuestra también que la circulación del campo vectorial a lo largo de la circunferencia se anula.
Para poder ilustrar esta fuerza se necesitan tomar soluciones de modelos de fluidos con una viscosidad positiva.

== º Curvas de nivel de la presión==
Observando las curvas de nivel de la presión podemos entender mejor las conclusiones a las que hemos ido llegando a lo largo del trabajo. En esta gráfica podemos observar los puntos de presión máxima <math>\theta=\frac{\pi}{2}</math> y <math>\theta=\frac{3\pi}{2}</math> y los puntos de presión mínima <math>\theta=0</math> y <math>\theta=\pi</math>.
Además, cabe destacar que, en las zonas en las que las curvas de nivel están mas juntas, hay mayor variación de la presión.

[[Archivo:Doce_mnj.jpg|800px|miniaturadeimagen|centro|Líneas de presión]]

Código Matlab:
{{matlab| codigo=

contour(x,y,P,100) % Curvas de nivel de la presión
axis([-5,5,-5,5]) % Región a dibujar

}}
[[Categoría:Teoría de Campos]]
[[Categoría:TC21/22]]