Diferencia entre revisiones de «Flujo de Poiseuille (GRUPO 10A)»

De MateWiki
Saltar a: navegación, buscar
(Líneas de corriente)
 
(No se muestran 156 ediciones intermedias de 2 usuarios)
Línea 1: Línea 1:
{{ TrabajoED |Flujo de Poiseuille (GRUPO 10A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC23/24|2023-24]] | Alba Prats Moreno <br/> Carlos Muñoz González <br/>Carla De Juan Merchán <br/> Rodrigo Prado Fornos <br/> Miguel Vela Gonçalves Cerejeira}}
+
{{ TrabajoED |Flujo de Poiseuille (GRUPO 10A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC23/24|2023-24]] | Alba Piedad Prats Moreno <br/> Carlos Muñoz González <br/>Carla De Juan Merchán <br/> Rodrigo Prado Fornos <br/> Miguel Vela Gonçalves Cerejeira}}
  
 
==Introducción==
 
==Introducción==
La ley de Poiseuille se utiliza para describir el flujo estacionario y laminar de un líquido incomprensible.  
+
La ley de Poiseuille se utiliza para describir el flujo estacionario y laminar de un líquido incompresible; es decir, la densidad del líquido solo cambia si se le aplica una presión determinada.
  
En el estudio de esta ley nos enfocamos en el flujo de un líquido incomprensible a través de una tubería cilíndrica cuyo radio es 2 centrada en el eje OZ.  
+
En el estudio de esta ley nos enfocamos en el flujo de un líquido incompresible a través de una tubería cilíndrica cuyo radio es 2 centrada en el eje OZ.  
  
La magnitud de este flujo viene determinada por el gradiente de presión y el radio de la propia tubería, para el desarrollo de este estudio hemos utilizado los programas Octave y Matlab y hemos trabajado en coordenadas cilíndricas.
+
La magnitud de este flujo viene determinada por el gradiente de presión y el radio de la propia tubería, teniendo en cuenta la viscosidad del líquido y la longitud de la tubería.  
  
 +
Para el desarrollo de este estudio hemos utilizado los programas Octave y Matlab y hemos trabajado en coordenadas cilíndricas.
  
  
==Mallado de la sección de la tubería==
 
Representación del mallado de dimensión 2, de la sección longitudinal del eje x=0. Consideraremos la región encerrada en las coordenadas (rho,z)=[0,3]x[0,10].
 
  
 +
[[Archivo: IntroG10Aa.png|300px|miniaturadeimagen|centro]]
  
[[Archivo: EJ1G10A.jpeg|400px|miniaturadeimagen|Mallado de la sección]]
+
==Mallado de la sección de la tubería==
 +
Representación del mallado de dimensión 2, de la sección longitudinal del eje x=0. Consideraremos la región encerrada en las coordenadas (ρ,z) = [0,3] x [0,10].
  
  
 +
[[Archivo: EJ1G10A.jpeg|300px|miniaturadeimagen|Mallado de la sección]]
 
  {{matlab|codigo=
 
  {{matlab|codigo=
 
  %1. Definimos los ejes
 
  %1. Definimos los ejes
Línea 31: Línea 33:
 
title ('Mallado de la sección');
 
title ('Mallado de la sección');
 
hold off
 
hold off
}}  
+
}}
  
 
==Ecuación de Navier-Stokes Estacionaria==
 
==Ecuación de Navier-Stokes Estacionaria==
La ecuación de Navier-Stokes describe como se mueve un fluido newtoniano. Esta herramienta es esencial para comprender como se comportan los fluidos en sistemas hidráulicos, como tuberías, canales...  
+
La ecuación de Navier-Stokes describe como se mueve un fluido newtoniano. El fluido newtoniano se caracteriza porque su deformación es proporcional al esfuerzo cortante (fuerza tangencial).
 +
 
 +
 
 +
Esta herramienta es esencial para comprender como se comportan los fluidos en sistemas hidráulicos, como tuberías, canales...  
 
Antes de sumergirnos en la demonstración es crucial establecer que el fluido es incompresible ya que esta ecuación rige el comportamiento de los fluidos newtonianos.
 
Antes de sumergirnos en la demonstración es crucial establecer que el fluido es incompresible ya que esta ecuación rige el comportamiento de los fluidos newtonianos.
En el enunciado nos proporcionan la siguiente igualdad:
+
 
<math>\overrightarrow{u}\left ( \rho ,\theta ,z \right )=f\left ( \rho  \right )\overrightarrow{ez}</math> y su presión <math>p\left ( x,y \right )=p1+\left ( p2-p1 \right )\left ( z-1 \right )</math>, donde <math>\mu</math> es el coeficiente de viscosidad, p1 es la presión en los puntos z = 1, p2 la presión en los puntos z = 3.
+
 
Para entender la formula es importante conocer a que se refiere cada termino, pues <math>\mu</math> es el coeficiente de viscosidad, p1 representa la presión para valores de z=1, y por consiguiente p2 lo hará para los valores de z=2.
+
Vamos a definir nuestra fuerza según su presión y superficie, relacionándolo con la viscosidad del medio. En el enunciado nos proporcionan la siguiente igualdad:
 +
 
 +
 
 +
<center><math>\overrightarrow{u}\left ( \rho ,\theta ,z \right )=f\left ( \rho  \right )\overrightarrow{ez}</math> </center>
 +
 
 +
y su presión:  <center><math>p\left ( x,y \right )=p1+\left ( p2-p1 \right )\left (\frac{z-1}{2 } \right )</math> </center>
 +
 
 +
 
 +
 
 +
Las variables que compreenden la ecuación de Navier-Stokes son las siguiente:
 +
 
 +
- <math>\mu</math> es el coeficiente de viscosidad.
 +
 
 +
- p1 representa la presión para valores de z=1.
 +
 
 +
- p2 lo hará para los valores de z=3.
 +
 
 +
 
 
El primer térmido de la ecuación será:
 
El primer térmido de la ecuación será:
<math>(\vec{u} \cdot \nabla) \vec{u})</math>
+
 
Si sustituimos por nuestro campo y lo desarrollamos veremos que nos queda el vector nulo por lo que el primer término de la ecuación de Navier-Stokes no interviene.
+
<center><math>(\vec{u} \cdot \nabla) \vec{u})</math></center>
 +
 
 +
 
 +
Si desarrollamos este primer término, como el fluido es incompresible, no hay fuerzas externas (fuentes ni sumideros) por lo que la divergencia va a ser nula. Consequentemente, podemos despreciar el primer término de la ecuación Navier-Stokes obteniendo:
 +
 
  
 
<center><math> \frac{1}{\rho } \frac{\partial }{\partial \rho } \left ( \rho \frac{\partial f(\rho )}{\partial \rho } \right ) = \frac{p_{2}-p_{1}}{\mu }</math></center>
 
<center><math> \frac{1}{\rho } \frac{\partial }{\partial \rho } \left ( \rho \frac{\partial f(\rho )}{\partial \rho } \right ) = \frac{p_{2}-p_{1}}{\mu }</math></center>
  
Esta última ecuación es una derivada segunda del campo requerido, por lo que debemos integrar dos veces, una vez múltiplicado por <math> \rho </math>
+
 
 +
 
 +
Esta ecuación representa la derivada segunda de f( <math> \rho </math>). Para obtener el valor de  <math> \rho </math> podemos hacer una integral doble y a continuación multiplicar por <math> \rho </math> ya que se trata de un gradiente en coordenadas cilíndricas:
 +
 
  
 
<center><math>f\left ( \rho  \right )=\frac{p2-p1}{4\mu}\rho ^{2} +Cln\left (\rho \right )+K</math>.</center>
 
<center><math>f\left ( \rho  \right )=\frac{p2-p1}{4\mu}\rho ^{2} +Cln\left (\rho \right )+K</math>.</center>
  
Para continuar debemos comprender el sistema, y saber que en <math> \rho=2 </math> la ecuación será 0 y para <math> \rho=2 </math> no tenderá a infinito, lo que nos deja lo siguiente:
 
  
<math>f\left (2\right )=0</math>,
+
Como se trata de una ecuación diferencial, vamos ahora a encontrar una solución particular para las constantes de integración.
 +
 
 +
Para encontrar la solución particular de la EDO, vamos a asumir el siguiente supuesto:
 +
 
 +
 
 +
La velocidad del fluido converge a 0 en las paredes de la tubería. Consequentemente, mayor será la velocidad del fluido cuando nos acercamos al centro de la tibería. Como el radio de la tubería es 2, el valor de la velocidad en los extremos es nula. 
 +
 
 +
 
 +
Es decir,  <center><math>f\left (2\right )=0</math>,
 +
<math>f\left ( 0 \right )=0</math> </center>,
 +
 
 +
 
 +
Vamos a proponer un sistemas de dos ecuaciones y obtener las constantes de integración:
 +
 
 +
 
 
<center><math>K=\frac{p1-p2}{\mu}</math> </center>
 
<center><math>K=\frac{p1-p2}{\mu}</math> </center>
 +
 +
<center><math>C=0</math> </center>
 +
 +
 +
Finalmente obtenemos la expresión f(<math> \rho </math>):
  
<math>f\left ( 0 \right )=0</math>,
 
<math>C=0</math>
 
  
 
<center><math>f\left ( \rho  \right )=\frac{1}{\mu }\left [ \frac{p2-p1}{4 }\rho ^{2} +p1-p2\right ]</math></center>
 
<center><math>f\left ( \rho  \right )=\frac{1}{\mu }\left [ \frac{p2-p1}{4 }\rho ^{2} +p1-p2\right ]</math></center>
  
Como ya se comentó al principio del apartado, el fluido es incompresible, para comprobarlo calcularemos la divergencia del campo, y al darnos 0 podemos concluir que la hipótesis inicial es valida, pues en caso de ser el resultado negativo significaría que se contrae y lo contrario en caso de resultado positivo.
+
 
 +
Para comprobar la condición de incomprensibilidad (el agua siempre ocupa el mismo volumen), vamos a comprobar que la divergencia del campo es nula.
 +
 
  
 
<center><math>\bigtriangledown \cdot \overrightarrow{u}=0=\frac{1}{\rho } \left\{ \frac{\partial }{\partial \rho }\left ( \rho u_{\rho } \right )+ \frac{\partial }{\partial \theta }(u_{\theta })+ \frac{\partial }{\partial z}(\rho u_{z})\right\}</math></center>
 
<center><math>\bigtriangledown \cdot \overrightarrow{u}=0=\frac{1}{\rho } \left\{ \frac{\partial }{\partial \rho }\left ( \rho u_{\rho } \right )+ \frac{\partial }{\partial \theta }(u_{\theta })+ \frac{\partial }{\partial z}(\rho u_{z})\right\}</math></center>
  
==Campo de presiones y campo de presiones==
 
Sustituyendo en la expresión del campo obtenemos:
 
  
<center><math>\vec{u}(\rho,\theta,z)=(\frac{-\rho ^{2}}{4 } +1)\vec{e_z}</math></center>
+
En conclusión, como f(p) no depende de "z", la divergencia es nula.
Para su representación se ha introducido en Matlab, el siguiente código:  
+
 
 +
==Campo de presiones y campo de velocidades==
 +
 
 +
 
 +
Vamos a dar como dato: <math> p1=4 </math>, <math> p2=1 </math> y <math> \mu=1 </math>.
 +
 
 +
 
 +
 
 +
===Campo de velocidades===
 +
 
 +
Anteriormente hemos obtenido la función de velocidad. Vamos a proceder a sustituir los valores dados:´
 +
 
 +
<center><math>\vec{u}(\rho,\theta,z)=(\frac{-3\rho ^{2}}{4} +3)\vec{e_z}</math></center>
 +
 
 +
Como podemos comprobar, a medida que nos acercamos al borde de la tubería, entra menos velocidad mientras que en el centro, hay mayor flujo de concentración de velocidades.
 +
 
 +
[[Archivo: Ejercicio3AG10A.jpeg|300px|miniaturadeimagen|Campo de velocidades.]]
 
{{matlab|codigo=
 
{{matlab|codigo=
x=0:0.2:2;
+
x=0:0.2:2; %Intervalo de 'ρ'
y=0:0.2:10;
+
z=0:0.2:10; %Intervalo de 'z'
[X,Y]=meshgrid(x,y);
+
[X,Z]=meshgrid(x,z);
ux=(-3./4).*xx.^2+3;
+
ux=(-3./4).*X.^2+3; %Dar valor a 'u'
uy=0.*yy;
+
uz=0.*zz;
 
hold on
 
hold on
quiver(xx,yy,ux,uy)
+
quiver(X,Z,ux,uz) %Representación
 
axis([0,5,0,12])
 
axis([0,5,0,12])
 
xlabel('ρ');
 
xlabel('ρ');
Línea 81: Línea 142:
 
hold off
 
hold off
 
view(2)
 
view(2)
title('CAMPO DE VELOCIDADES')
+
title('Campo de velocidades')
 
}}
 
}}
  
 +
===Campo de presiones===
  
{{matlab|codigo=
+
La expresión del campo de presiones viene dada como:
x=0:0.2:2;
+
z=0:0.2:10;
+
[X,Z]=meshgrid(x,z);
+
figure(1);
+
p=4-3/2*(Z-1);
+
surf(X,Z,p);
+
colorbar;
+
view(2);
+
axis([0,3,0,10]);
+
title('CAMPO DE PRESIONES');
+
xlabel('ρ');
+
ylabel('z');
+
}}
+
  
==Velocidad máxima del fluido y módulo de velocidad==
+
<center><math> p(x,y,z)=p_1+\frac{(p_2-p_1)(z-1)}{2}</math></center>
  
Los puntos de velocidad máxima se obtienen igualando a 0, la primera derivada parcial del campo:
+
Sustituyendo los valores:  
  
<center><math>\vec{u}(\rho,\theta,z)=(\frac{-\rho ^{2}}{4 } +1)\vec{e_z}</math></center>
+
<center><math>p_1=4 \ {,} \ p_2=1 \ {,}  \ μ=1</math> </center>
  
donde:
+
Obtenemos por tanto:
 +
<center><math> p(x,y)=4-\frac{3z+3}{2}=\frac{11}{2}-\frac{3z}{2}</math></center>
  
<center><math> \frac{\partial \vec{u}} {\partial \rho}=\frac{-\rho}{2}\vec{e_z} </math></center>
 
<center><math> \frac{\partial \vec{u} }{\partial \theta}=0 </math></center>
 
<center><math> \frac{\partial \vec{u} }{\partial z}=0 </math></center>
 
  
Si igualamos a 0, <math> \frac{-\rho}{2}=0</math> obtenemos que el punto donde se produce la máxima velocidad es <math> \rho =0 </math>, como hemos podido observar en el gráfico anterior del campo de velocidades.
+
Como podemos comprobar, la presión no aumenta con el radio, sino con el parámetro "z".
  
A continuación, se muestra la curva de la variación de la velocidad del fluido en función del radio de la sección del tubo circular. Como se puede apreciar en el gráfico siguiente a medida que aumenta el radio, disminuye la velocidad, tal y como hemos reflejado en apartados anteriores.
+
[[Archivo:CdPG10Aa.jpg|250px|miniaturadeimagen|Campo de presiones]]
  
[[Archivo:modvel.png|400px|miniaturadeimagen| Módulo de la velocidad]]
 
  
{{matlab|codigo=
 
rho=0:0.1:2;
 
f=abs((-(rho.^2/4)+1);
 
plot(rho,f);
 
title('Módulo del campo de velocidades');
 
axis([0,2,0,3.5])
 
}}
 
  
==CAMPO DE PRESIONES==
 
El campo de presiones al que está sometido nuestro fluido es un campo escalar, con las siguiente expresión:
 
<center><math> p(x,y,z)=p_1+(p_2-p_1)(z-1)</math></center>
 
Para representarlo gráficamente se consideran nuevamente los valores:
 
<math>p_1=2 \ {,} \  p_2=1 \ {,}  \  μ=1</math>
 
.Obteniendo, por tanto:
 
<center><math> p(x,y)=2-(z-1)=3-z</math></center>
 
  
Implementado el código siguiente en Octave:
 
[[Archivo:cdp27.png|400px|miniaturadeimagen|Campo de presiones]]
 
 
{{matlab|codigo=
 
{{matlab|codigo=
y=0:0.05:2;
+
z=0:0.1:10; %Altura 'z'
z=0:0.05:10;
+
f=(-3)*z/2+(11/2); %Funcion del campo de presiones
[Y,Z]=meshgrid(y,z);  
+
plot(z,f,'r')
figure (1)
+
xlabel('Altura(z)');
p=3-Z;
+
ylabel('Presión(p)');
surf(Y,Z,p)
+
title(' Campo de presiones');
view(2)
+
axis([0,5,0,12])
+
 
}}
 
}}
  
 +
== Líneas de corriente ==
  
 +
Para dibujar las líneas de corriente del campo <math>\overrightarrow{u}</math>, debemos tener en cuenta que estas son tangenes a <math>\overrightarrow{u}</math> en cada apunto.
  
 +
Procedemos a calcular <math>\overrightarrow{v}</math> que es ortogonal a <math>\overrightarrow{u}</math> ya que se comprueba que:
 +
<math>\overrightarrow{v}=\overrightarrow{e_{\theta }}\cdot \overrightarrow{u}</math>.
  
  
 +
Como hemos hallado anteriormente, la divergencia de <math>\overrightarrow{u}</math> es nula. Consequentemente, el rotacional del campo <math>\overrightarrow{v}</math> es nulo también.
  
  
 
+
La función de corriente o potencial escalar viene definido como:  <math>\psi </math>,(<math>\bigtriangledown \psi =\overrightarrow{v}</math>)  
 
+
 
+
 
+
 
+
 
+
<center> <font color="D2 69 1E">'''INTERPRETACIÓN'''</font> </center>
+
 
+
Como se puede apreciar en el gráfico a medida que el fluido avanza por el canal la presión va disminuyendo. Las presiones más altas están representadas con colores cálidos, y las presiones más bajas con colores fríos. Esto se debe a que para mantener el caudal de un fluido viscoso estable debe mantenerse una diferencia de presiones entre las paredes del canal, esta diferencia de presión es necesaria debida a la fuerza de arrastre o frenada que ejerce el canal sobre la capa de fluido en contacto con él. Estas fuerzas de arrastre o de frenado se denominan fuerzas viscosas. El resultado, hace que la velocidad del fluido no sea constante a lo largo del canal siendo mayor cerca del centro y menor cerca de las paredes.
+
=== LÍNEAS DE CORRIENTE ===
+
Vamos a dibujar las líneas de corriente del campo <math>\overrightarrow{u}</math>, es decir, las líneas que son tangenes a <math>\overrightarrow{u}</math> en cada apunto, es necesario calcular el campo <math>\overrightarrow{v}</math> que en cada punto es ortogonal a <math>\overrightarrow{u}</math>, Siendo
+
<math>\overrightarrow{v}=\overrightarrow{e_{\theta }}\cdot \overrightarrow{u}</math>.
+
El campo <math>\overrightarrow{v}</math> es irrotacional por ser nula la divergencia de <math>\overrightarrow{u}</math> y ,además, tiene un potencial escalar <math>\psi </math>,(<math>\bigtriangledown \psi =\overrightarrow{v}</math>), que se conoce como función de corriende de <math>\overrightarrow{u}</math>.
+
  
 
<center><math>\overrightarrow{v}=\begin{vmatrix}
 
<center><math>\overrightarrow{v}=\begin{vmatrix}
Línea 173: Línea 194:
 
\end{vmatrix}=f\left ( \rho  \right )\overrightarrow{e_{\rho }} \rightarrow \overrightarrow{v}=f\left ( \rho  \right )\overrightarrow{e_{\rho }}</math>.</center>
 
\end{vmatrix}=f\left ( \rho  \right )\overrightarrow{e_{\rho }} \rightarrow \overrightarrow{v}=f\left ( \rho  \right )\overrightarrow{e_{\rho }}</math>.</center>
  
Sustituyendo <math>f\left ( \rho  \right )</math> obtenemos:
+
 
 +
Sustituyendo <math>f\left ( \rho  \right )</math>:
  
 
<center><math>\overrightarrow{v}=\left ( \frac{p2-p1}{4\mu }\rho ^{2} +\frac{p1-p2}{\mu }\right )\overrightarrow{e_{\rho }}</math>.</center>
 
<center><math>\overrightarrow{v}=\left ( \frac{p2-p1}{4\mu }\rho ^{2} +\frac{p1-p2}{\mu }\right )\overrightarrow{e_{\rho }}</math>.</center>
  
Para calcular <math>\psi </math>, obtenemos su gradiente e igualamos al campo <math>\overrightarrow{v}</math>, <math>\bigtriangledown \psi =\overrightarrow{v}</math>.
 
  
 
<center><math>\bigtriangledown \psi=\frac{d\psi }{d\rho }\overrightarrow{e_{\rho }}+\frac{1}{\rho }\frac{d\psi }{d\theta }\overrightarrow{e_{\theta }}+\frac{d\psi }{dz}\overrightarrow{e_{z}}=\overrightarrow{v}</math>.</center>
 
<center><math>\bigtriangledown \psi=\frac{d\psi }{d\rho }\overrightarrow{e_{\rho }}+\frac{1}{\rho }\frac{d\psi }{d\theta }\overrightarrow{e_{\theta }}+\frac{d\psi }{dz}\overrightarrow{e_{z}}=\overrightarrow{v}</math>.</center>
  
Por ende: <math>\frac{d\psi }{d\rho }=\frac{p2-p1}{4\mu }\rho ^{2} +\frac{p1-p2}{\mu }\ </math>.
 
  
Tras integrar nos queda finalmente que <math>\psi=\frac{p2-p1}{24\mu }\rho ^{2} +\frac{p1-p2}{\mu }\rho </math>.
 
  
Sustituyendo en la anterior función los valores asignados, el resultado de la función potencial será: <math> \psi= -\frac{1}{24} \rho ^{2} + \rho</math>.
+
Obtenemos la relación: <center><math>\frac{d\psi }{d\rho }=\frac{p2-p1}{4\mu }\rho ^{2} +\frac{p1-p2}{\mu }\ </math>.</center>
 +
 
 +
 
 +
Integrando el gradiente del potencial escalar:    <center><math>\psi=\frac{p2-p1}{12\mu }\rho ^{3} +\frac{p1-p2}{\mu }\rho </math>.</center>
 +
 
 +
 
 +
 
 +
El resultado de la función potencial quedará: <center><math> \psi= -\frac{1}{4} \rho ^{3} + 3\rho</math>.</center>
 +
 
 +
 
 +
[[Archivo:LineasG10Aa.jpg|250px|miniaturadeimagen|Líneas de corriente]]
  
Con ayuda del programa Matlab se muestran las líneas de corriente de la siguiente manera:
 
[[Archivo:lineascorriente27a.png|400px|miniaturadeimagen|Líneas de corriente]]
 
 
{{matlab|codigo=
 
{{matlab|codigo=
x=0:0.2:2;
+
x=0:0.2:2; %Definimos 'ρ'
y=0:0.2:10;
+
z=0:0.2:10; %Definimos 'z'
[X,Y]=meshgrid(x,y);
+
[X,Z]=meshgrid(x,z);
lineas=(-1./24).*X.^2+X;
+
lineas=(-1./4).*X.^3+3.*X; %Definimos campo escalar 
contour(X,Y,lineas);
+
contour(X,Z,lineas);  
 
axis([0,3,0,10]);
 
axis([0,3,0,10]);
 
colorbar
 
colorbar
Línea 200: Línea 227:
 
}}
 
}}
  
 +
==Velocidad máxima del fluido y módulo de velocidad==
  
 +
Para obtener los máximos y mínimos de la velocidad vamos a optimizar la función de velocidad. Para ello, derivamos la misma y la igualamos a 0:
  
  
 +
<center><math>\vec{u}(\rho,\theta,z)=(\frac{-3\rho ^{2}}{4 } +3)\vec{e_z}</math></center>
  
  
 +
donde:
  
 +
<center><math> \frac{\partial \vec{u}} {\partial \rho}=\frac{-3\rho}{2}\vec{e_z}</math></center>
 +
<center><math> \frac{\partial \vec{u} }{\partial \theta}=0 </math></center>
 +
<center><math> \frac{\partial \vec{u} }{\partial z}=0 </math></center>
  
  
 +
Consequentemente, el punto donde la velocidad es máxima es en <math> \rho =0 </math>
  
 +
Como podemos observar en la gráfica la velocidad de nuestro fluido no es constante a lo largo de la tubería. De hecho, disminuye a medida que nos acercamos a los bordes de la tubería. Esto se debe a las fuerzas viscosas potenciadas por la diferencia de presión que debe mantener el fluido para ser viscoso.
 +
 +
 +
[[Archivo: Ejercicio3BG10Aa.jpg|270px|miniaturadeimagen|Módulo velocidad de fluido]]
  
<center> <font color="D2 69 1E">'''INTERPRETACIÓN'''</font> </center>
 
Como se puede observar en el gráfico anterior, estas líneas son tangentes al campo vectorial <math>\vec{u}</math> en cada punto, quedando demostrado que son las líneas de corriente de <math>\vec{u}</math>.
 
  
==ROTACIONAL==
 
Para la realización del cálculo del rotacional del campo, utilizaremos la siguiente fórmula:
 
  
<center> <big> <math>\nabla\times\vec u= \frac{1}{\rho}\begin{vmatrix} \vec{e_\rho} & \rho \vec{e_\theta} & \vec{e_z} \\ \frac{ \partial}{\partial \rho} & \frac{\partial }{\partial \theta} & \frac{\partial }{\partial z}\\ u_\rho & \rho u_\theta & u_z \end{vmatrix}</math></big></center>
 
  
Sustituyendo en la fórmula:
 
<center> <big> <math>\nabla\times\vec u= \frac{1}{\rho}\begin{vmatrix} \vec{e_\rho} & \rho \vec{e_\theta} & \vec{e_z} \\ \frac{ \partial}{\partial \rho} & \frac{\partial }{\partial \theta} & \frac{\partial }{\partial z}\\ 0 & 0 & (\frac{-\rho ^{2}}{4 } +1)\frac{\cdot(p_1-p_2)}{μ} \end{vmatrix}</math></big></center>
 
Obtenemos que <math>\nabla\times\vec{u}=(\frac{-3\rho}{4}+\frac{1}{\rho})\frac{(p_1-p_2)}{μ}\vec{e_\theta}</math>, y dando los siguientes valores:
 
<center><math>p_1=2 \ {,} \  p_2=1 \ {,}  \  μ=1</math></center>
 
Tenemos que <math>\nabla\times\vec u=(\frac{-3\rho}{2}+\frac{2}{\rho})\vec{e_\theta}</math>, siendo el módulo del rotacional,<math>\left | \nabla\times\vec u \right |= (\frac{-3\rho}{2}+\frac{2}{\rho}) </math>, como se observa en la expresión,  depende del valor de <math>\rho</math>. Tomando los valores mínimos en los extremos del intervalo en el que está definido <math>\rho</math>, <math> [0,2]</math>.
 
  
Para visualizar esto de manera gráfica, nuevamente a Octave,
 
[[Archivo:Rotacional27a.png|400px|miniaturadeimagen|Rotacional del campo velocidad]]
 
 
{{matlab|codigo=
 
{{matlab|codigo=
x=0:0.1:2;
+
x=0:0.05:2; %Definimos rho
y=0:0.1:10;
+
f=abs((-3/4)*x.^2+3); %Definimos funcion de velocidad
[xx,yy]=meshgrid(x,y);
+
plot(x,f,'r');
rot=abs(((-y.^2)./4)+1);
+
title('Módulo de la velocidad')
surf(xx,yy,rot)
+
xlabel('Radio')
colorbar
+
ylabel('Velocidad')
view(2)
+
axis([0,3,0,5])
axis([0,5,0,12])
+
 
}}
 
}}
  
 +
==Rotacional==
 +
Para hallar el rotacional, vamos a proceder con el cálculo:
  
 +
<center> <big> <math>\nabla\times\vec u= \frac{1}{\rho}\begin{vmatrix} \vec{e_\rho} & \rho \vec{e_\theta} & \vec{e_z} \\ \frac{ \partial}{\partial \rho} & \frac{\partial }{\partial \theta} & \frac{\partial }{\partial z}\\ u_\rho & \rho u_\theta & u_z \end{vmatrix}</math></big></center>
  
 +
Sustituyendo en la fórmula:
 +
<center> <big> <math>\nabla\times\vec u= \frac{1}{\rho}\begin{vmatrix} \vec{e_\rho} & \rho \vec{e_\theta} & \vec{e_z} \\ \frac{ \partial}{\partial \rho} & \frac{\partial }{\partial \theta} & \frac{\partial }{\partial z}\\ 0 & 0 & (\frac{-3\rho ^{2}}{4 } +3)\end{vmatrix}</math></big></center>
  
 +
Obtenemos que <math>\nabla\times\vec{u}=(\frac{3\rho}{2})\vec{e_\theta}</math>,
  
  
 +
Para visualizar esto de manera gráfica, nuevamente a Octave,
 +
[[Archivo:RotG10Aa.jpg|250px|miniaturadeimagen|Rotacional del campo velocidad]]
  
 +
{{matlab|codigo=
 +
x=0:0.1:2; %Definimos 'rho'
 +
z=0:0.1:10; %Definimos 'z'
 +
[X,Z]=meshgrid(x,z);
 +
rot=abs(3.*X./2); %Calculamos el rotacional
 +
surf(X,Z,rot)
 +
colorbar
 +
view(2)
 +
axis([0,5,0,12])
 +
}}
  
 +
==Campo de temperaturas==
  
 +
En el enunciado nos viene dado la función Temperatura en coordenadas cilíndricas:
  
 +
<center><math> T(ρ{,}θ{,}z)=log(1+ρ)e^{-(z-2)^2}</math></center>
  
<center> <font color="D2 69 1E">'''INTERPRETACIÓN'''</font> </center>
 
El rotacional muestra la tendencia del campo a inducir rotación alrededor de un punto. Por otra parte, el rotacional caracteriza la rotación de un fluido, por lo que en los extremos del canal el efecto de giro será máximo, particularizando en <math>\rho=2</math> tendrá el valor máximo en el sentido positivo de <math>\vec{e_z}</math>.
 
En el gráfico los puntos con mayor tendencia a la rotación aparecen representados con colores cálidos, amarillos hacia verde azulado, y los puntos con menor tendencia a la rotación aparecen representados con colores fríos, azules. Como se puede apreciar en el gráfico aparece la tendencia del rotacional en torno a las paredes del canal <math>\rho=2</math>.
 
 
== CAMPO DE TEMPERATURAS==
 
  
La temperatura a la que está sometido el fluido viene determinada por el campo escalar:
+
En la gráfica podemos observar como la temperatura es máxima para p=2, z=2; ya que los valores p y z son proporcionales a la temperatura.
<center><math> T(ρ{,}θ{,}z)=1+(ρ-1/2)^2 e^{-(z-1)^2}</math></center>
+
  
Como se observa, la expresión del campo está en coordenadas porlares. Teniendo en cuenta que las relación para pasar de coordenadas polares a cartesianas es:
 
: <math>x=ρcosθ</math>
 
: <math>y=ρsenθ</math>
 
  
Pero como nosotros estamos en los ejes <math>\{\vec{j} {,} \vec{k} \}</math>, usaremos:
+
[[Archivo:CdTG10Aa.png|320px|miniaturadeimagen| Campo de temperaturas]]
 
+
{{Matlab|codigo=
: <math>y=ρcosθ</math>
+
p=0:0.01:2; %Definimos p
: <math>z=ρsenθ</math>
+
z=-2:0.05:10; %Definimos z
Podemos expresar el campo en coordenadas cartesianas, que es el sistema en el cual estamos trabajando, resultando:
+
[P,Z]=meshgrid(p,z);
<center><math> T(y{,}z)=1+(y^2+z^2-1/2)^2 e^{-(z-1)^2}</math></center>
+
figure (1)  
 
+
a=log(1+P)*exp.^-(-Z-2).^2 %Funcion de temperatura
 
+
pcolor(Y,Z,p);
<font color="80 80 00">''' REPRESENTACIÓN DEL CAMPO DE TEMPERATURAS '''</font>
+
shading flat
 
+
El campo de temperaturas es un campo escalar que representa la distribución espacial de la temperatura asociando un valor a cada punto del espacio. Depende de la posición del punto, y del tiempo (t). En nuestro caso, el campo de temperaturas como ya se a expuesto viene dado por la expresión  <math> T(y{,}z)=1+(y^2+z^2-1/2) e^{-((z-1)^2)}</math>, por lo tanto no depende del tiempo, se dice que es estacionario, sino exclusivamente de las componentes espaciales <math>(y{,}z)</math>
+
 
+
Para la representación se ha implementado el siguiente código en Octave:
+
[[Archivo:Temperatura27-1.PNG|400px|miniaturadeimagen| Campo de temperaturas del fluido]]
+
{{matlab|codigo=
+
y=0:0.05:8;
+
z=-2:0.05:10;
+
[Y,Z]=meshgrid(y,z);
+
figure (1)  
+
p=1+(((Y.^2+Z.^2-(1/2)).^2).*exp(-((Z-1).^2)));
+
pcolor(Y,Z,p);
+
shading flat
+
grid on
+
axis([0,8,-1,10]);
+
colorbar
+
figure(2)
+
contour(Y,Z,p,10,'k');
+
 
grid on
 
grid on
 
axis([0,8,-1,10]);
 
axis([0,8,-1,10]);
 +
10 colorbar
 +
11 figure(2)
 +
12 contour(Y,Z,a,10,'k');
 +
13 grid on
 +
14 axis([0,8,-1,10]);
 
}}
 
}}
  
  
 +
===Curvas de nivel===
  
  
<center><font color="D2 69 1E">'''INTERPRETACIÓN'''</font> </center>
+
Las curvas de nivel son líneas que concentran puntos con la misma temperatura. Se comprueba que tienen forma simétrica y logarítmica.
En el gráfico se puede observar que la temperatura es mayor cuándo se representa con colores cálidos, esto sucede en los valores próximos a <math>z=1</math>.
+
  
<font color="80 80 00">''' REPRESENTACIÓN DE LAS CURVAS DE NIVEL '''</font>
 
[[Archivo:Temperatura27-2.PNG|400px|miniaturadeimagen|centro| Curvas de nivel de campo]]
 
  
El código utilizado para la representación de este gráfico se incluyó con el de la representación del campo de temperaturas.
 
  
<center><font color="D2 69 1E">'''INTERPRETACIÓN'''</font> </center>
+
[[Archivo:CdNG10Aa.jpg|300px|miniaturadeimagen| Curvas de nivel de campo de temperatura]]
 +
{{Matlab|codigo=
 +
x=0:0.05:2; %Definimos 'rho'  
 +
z=0:0.05:1; %Definimos 'z'
 +
[X,Z]=meshgrid(x,z);
 +
figure (1)
 +
p=log(1+P)*exp.^-(-Z-2).^2 %Funcion de temperatura
 +
[pX,pZ]=gradient(p); %Funcion gradiente
 +
hold on
 +
quiver(X,Z,pX,pZ)
 +
axis([0,8,-1,2]);
 +
shading flat
 +
grid on
 +
hold off
 +
}}
  
En el gráfico aparecen representadas las curvas de nivel, las curvas de nivel varían de manera geométrica estando más próximas entre sí cuando la temperatura aumenta, en torno a los puntos <math>z=1</math> y más alejadas cuando la temperatura disminuya.
+
==Gradiente de temperatura==
  
La representaciónes anteriores permite determinar que en el punto <math>(y,z)= (8,1)</math> la temperatura es máxima.
+
Teniendo la función de temperatura, calculamos su gradiente:
  
===GRADIENTE DEL CAMPO DE TEMPERATURAS===
+
<center><math>\nabla T(ρ,θ,z)=(\frac{e^{-(z-2)^2}}{1+ρ)})\vec{e_\rho}-(2log(ρ+1)(z-2){e^{-(z-2)^2}})\vec{e_z}</math></center>
El análisis de la variación de la temperatura a lo largo del canal se realiza analizando el gradiente <math>\nabla T</math>, en nuestro caso:
+
  
<center><math>\nabla T(y,z)=(4y^3e^{-(z-1)^2}+4yz^2e^{-(z-1)^2}-2ye^{-(z-1)^2})\vec{j}</math> <math>+(4z^3e^{-(z-1)^2}-2z^4e^{-(z-1)^2}+4zy^2e^{-(z-1)^2}-4y^2z^2e^{-(z-1)^2} + 2ze^{-(z-1)^2}-2ze^{-(z-1)^2}) \vec{k}</math></center>
 
  
Como observación, el gradiente analíticamente se calcula como <math>\frac{\partial T(y,z)}{\partial y}</math> siendo ésta la compontente <math>\vec{j}</math> del campo vectorial gradiente y <math>\frac{\partial T(y,z)}{\partial z }</math> la componente <math>\vec{k}</math>.
+
La gráfica del gradiente de la temperatura, se ha representado utilizando el siguiente código:
  
<font color="80 80 00">''' REPRESENTACIÓN GRÁFICA DEL GRADIENTE'''</font>'''
 
  
El gráfico del gradiente de la temperatura, se ha representado utilizando el siguiente código:
+
[[Archivo:GRadienteG10Aa.jpg|300px|miniaturadeimagen|Gradiente de la temperatura|derecha]]
[[Archivo:Temperatura27-4.PNG|400px|miniaturadeimagen| Campo vectorial gradiente]]
+
 
 
{{Matlab|codigo=
 
{{Matlab|codigo=
y=0:0.05:8;
+
rho=0:0.1:2; %Definimos 'rho'
z=0:0.05:1;
+
z=0:0.1:10; %Definimos 'z'
[Y,Z]=meshgrid(y,z);
+
[RHO,Z]=meshgrid(rho,z);
figure (1)
+
figure(1)
p=1+((Z.^2).*exp(-((Y.^2)+(Z.^2)-1/2).^2));
+
T=log(1+RHO).*exp(-(Z-2).^2); %Funcion temperatura
[pY,pZ]=gradient(p);
+
[TRHO,TZ]=gradient(T); %Funcion gradiente
 
hold on
 
hold on
quiver(Y,Z,pY,pZ)
+
quiver(RHO,Z,TRHO,TZ)
axis([0,8,-1,2]);
+
xlabel('rho');
 +
ylabel('z');
 +
axis([0,2,0,10]);
 +
title('Gradiente de Temperatura')
 
shading flat
 
shading flat
 
grid on
 
grid on
Línea 334: Línea 370:
 
}}
 
}}
  
Código Octave utilizado:
+
[[Archivo:GRayCurvG10Aa.jpg|300px|miniaturadeimagen|Gradiente y curvas de nivel|derecha]]
[[Archivo:Temperatura27-5.PNG|400px|thumb|Gradiente ortogonal a las curvas de nivel]]
+
 
 
{{Matlab|codigo=
 
{{Matlab|codigo=
y=0:0.05:8;
+
rho=0:0.1:2; %Definimos 'rho'
z=0:0.05:10;
+
z=0:0.1:10; %Definimos 'z'
[Y,Z]=meshgrid(y,z);
+
[RHO,Z]=meshgrid(rho,z);
figure (1)
+
figure(1)
p=1+(((Y.^2+Z.^2-(1/2)).^2).*exp(-((Z-1).^2)));
+
T=log(1+RHO).*exp(-(Z-2).^2); %Definimos la funcion de temperatura
[pY,pZ]=gradient(p);
+
[TRHO,TZ]=gradient(T); %Definimos la funcion gradiente
 
hold on
 
hold on
quiver(Y,Z,pY,pZ)
+
quiver(RHO,Z,TRHO,TZ)
contour(Y,Z,p,'k');
+
contour(RHO,Z,T,'k') %Colocamos las curvas de nivel superpuestas
axis([0,8,-1,10]);
+
xlabel('rho');
 +
ylabel('z');
 +
axis([0,2,0,10]);
 +
title('Gradiente de Temperatura y Curvas de nivel')
 
shading flat
 
shading flat
 
grid on
 
grid on
Línea 355: Línea 394:
  
  
 +
Como podemos verificar, las curvas de nivel son ortogonales al gradiente de la temperatura. Esto se debe a que el gradiente es el responsable de la dirección de crecimiento y decrecimiento del campo.
  
 +
Por otro lado, podemos comprobar que si el gradiente es mayor, la temperatura tendrá una variación más rápida.
  
 +
==Caudal==
  
 +
El caudal es la cantidad de un fluido que atraviesa la tubería en un determinado tiempo. Para ello debemos tener en cuenta el flujo de volumen por unidad de tiempo. Vamos a representarlo de la siguiente forma:
  
 +
<center><math>Q=\int_{S}\vec{v} \cdot d\vec{S}</math></center>
  
 +
donde <math>\vec{v}</math> es el campo de velocidades.
  
 +
El vector normal será perpendicular a la superficie por lo que hemos visto anteriormente de las curvas y gradientes.
  
 +
El campo de velocidades finalmente será:
  
Con este gráfico afirmamos que las curvas de nivel son ortogonales al gradiente de la temperatura.
+
<center><math>\vec{u}(\rho,\theta,z)=(-\frac{3}{4}\rho ^{2} +3)\overrightarrow{e_{z}}</math></center>,  
 
+
<center><font color="D2 69 1E">'''INTERPRETACIÓN'''</font> </center>
+
 
+
Para estudiar la variación de temperatura a lo largo del canal analizamos el gradiente que indica la direccion de crecimiento de la función temperatura en cada uno de los puntos. Como ya se ha mencionado antes, junto al gradiente aparecen las curvas de nivel de la temperatura. Se puede comprobar que en torno al punto <math>z=1</math> el módulo del gradiente aumenta, lo que indica que la temperatura tendrá una variación más rápida en esta zona. Se aprecia claramente como el gradiente es ortogonal a las curvas de nivel, esto se debe a que el gradiente  indica la dirección de crecimiento y sentido de los valores del campo en cada punto y las curvas de nivel son el lugar geométrico de los puntos equipotenciales.
+
 
+
==CAUDAL==
+
El caudal es el volumen de un fluido que pasa a través de la tubería estudiada por unidad de tiempo y se calcula mediante la integral de superficie siguiente:
+
 
+
<center><math>\int_{S}\vec{v} \cdot d\vec{S}</math></center>
+
 
+
donde <math>\vec{v}</math> es el campo de velocidades. Al tratarse de una integral de superficie hay que tener en cuenta el vector normal que es perpendicular a la superficie.
+
 
+
En nuestro caso, el campo de velocidades es:
+
 
+
<center><math>\vec{u}(\rho,\theta,z)=(-\frac{1}{4}\frac{\rho ^{2}}{\mu }+ 1)\overrightarrow{e_{z}}</math></center>,  
+
 
+
Esa expresión resulta de sustituir los valores siguientes en la expresión incial del campo:
+
 
+
<center><math>p_1=2 \ {,} \  p_2=1 \ {,}  \  μ=1</math></center>
+
 
+
La profundidad del caudal es de 1 metro, por lo que tras parametrizar realizamos el calculo del caudal.
+
  
'''CÁLCULO DEL CAUDAL'''
+
Finalmente resolviendo la integral doble obtendremos:
  
<math>Q\left ( m/s \right )=\int_{S}^{}\overrightarrow{u}\overrightarrow{n}dS=\int_{S}^{}\left ( \frac{1}{4}\frac{\rho ^{2}}{\mu }+\frac{\rho}{\mu }  \right )\overrightarrow{e_{z}}\cdot \overrightarrow{e_{z}}
+
<math>Q\left ( m/s \right )=\int_{S}^{}\overrightarrow{u}\overrightarrow{n}dS=\int_{S}^{}\left (\frac{-3\rho^{2}}{4}+{3} \right )\overrightarrow{e_{z}}\cdot \overrightarrow{e_{z}}
dS=\int_{0}^{2\Pi }\int_{0}^{2}\left ( \frac{1}{4}\frac{\rho ^{2}}{\mu }+\frac{\rho}{\mu }  \right )d_{\rho }d_{\theta }=11.52\left ( m/s \right )</math>.
+
dS=\int_{0}^{2\Pi }\int_{0}^{2}\left ( \frac{-3\rho^{2}}{4}+{3}  \right )d_{\rho }d_{\theta }=8{\Pi}=25.133\left ( m/s \right )</math>.
  
 
[[Categoría:Teoría de Campos]]
 
[[Categoría:Teoría de Campos]]
[[Categoría:TC22/23]]
+
[[Categoría:TC23/24]]

Revisión actual del 23:14 9 dic 2024

Trabajo realizado por estudiantes
Título Flujo de Poiseuille (GRUPO 10A)
Asignatura Teoría de Campos
Curso 2023-24
Autores Alba Piedad Prats Moreno
Carlos Muñoz González
Carla De Juan Merchán
Rodrigo Prado Fornos
Miguel Vela Gonçalves Cerejeira
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


1 Introducción

La ley de Poiseuille se utiliza para describir el flujo estacionario y laminar de un líquido incompresible; es decir, la densidad del líquido solo cambia si se le aplica una presión determinada.

En el estudio de esta ley nos enfocamos en el flujo de un líquido incompresible a través de una tubería cilíndrica cuyo radio es 2 centrada en el eje OZ.

La magnitud de este flujo viene determinada por el gradiente de presión y el radio de la propia tubería, teniendo en cuenta la viscosidad del líquido y la longitud de la tubería.

Para el desarrollo de este estudio hemos utilizado los programas Octave y Matlab y hemos trabajado en coordenadas cilíndricas.


centro

2 Mallado de la sección de la tubería

Representación del mallado de dimensión 2, de la sección longitudinal del eje x=0. Consideraremos la región encerrada en las coordenadas (ρ,z) = [0,3] x [0,10].


Mallado de la sección
%1. Definimos los ejes
rho=0:0.2:2; 
z=0:0.2:10; 
%2. Definimos mallado en 2 dimensiones
[xx,yy]=meshgrid(rho,z); %Mallado XY.
hold on
mesh(xx,yy,0.*xx); %Representamos la tubería
axis([0,3,0,10]); %Rango de los ejes
xlabel('ρ') ;
ylabel('z') ;
view(2);
title ('Mallado de la sección');
hold off


3 Ecuación de Navier-Stokes Estacionaria

La ecuación de Navier-Stokes describe como se mueve un fluido newtoniano. El fluido newtoniano se caracteriza porque su deformación es proporcional al esfuerzo cortante (fuerza tangencial).


Esta herramienta es esencial para comprender como se comportan los fluidos en sistemas hidráulicos, como tuberías, canales... Antes de sumergirnos en la demonstración es crucial establecer que el fluido es incompresible ya que esta ecuación rige el comportamiento de los fluidos newtonianos.


Vamos a definir nuestra fuerza según su presión y superficie, relacionándolo con la viscosidad del medio. En el enunciado nos proporcionan la siguiente igualdad:


[math]\overrightarrow{u}\left ( \rho ,\theta ,z \right )=f\left ( \rho \right )\overrightarrow{ez}[/math]
y su presión:
[math]p\left ( x,y \right )=p1+\left ( p2-p1 \right )\left (\frac{z-1}{2 } \right )[/math]


Las variables que compreenden la ecuación de Navier-Stokes son las siguiente:

- [math]\mu[/math] es el coeficiente de viscosidad.

- p1 representa la presión para valores de z=1.

- p2 lo hará para los valores de z=3.


El primer térmido de la ecuación será:

[math](\vec{u} \cdot \nabla) \vec{u})[/math]


Si desarrollamos este primer término, como el fluido es incompresible, no hay fuerzas externas (fuentes ni sumideros) por lo que la divergencia va a ser nula. Consequentemente, podemos despreciar el primer término de la ecuación Navier-Stokes obteniendo:


[math] \frac{1}{\rho } \frac{\partial }{\partial \rho } \left ( \rho \frac{\partial f(\rho )}{\partial \rho } \right ) = \frac{p_{2}-p_{1}}{\mu }[/math]


Esta ecuación representa la derivada segunda de f( [math] \rho [/math]). Para obtener el valor de [math] \rho [/math] podemos hacer una integral doble y a continuación multiplicar por [math] \rho [/math] ya que se trata de un gradiente en coordenadas cilíndricas:


[math]f\left ( \rho \right )=\frac{p2-p1}{4\mu}\rho ^{2} +Cln\left (\rho \right )+K[/math].


Como se trata de una ecuación diferencial, vamos ahora a encontrar una solución particular para las constantes de integración.

Para encontrar la solución particular de la EDO, vamos a asumir el siguiente supuesto:


La velocidad del fluido converge a 0 en las paredes de la tubería. Consequentemente, mayor será la velocidad del fluido cuando nos acercamos al centro de la tibería. Como el radio de la tubería es 2, el valor de la velocidad en los extremos es nula.


Es decir,
[math]f\left (2\right )=0[/math], [math]f\left ( 0 \right )=0[/math]
,


Vamos a proponer un sistemas de dos ecuaciones y obtener las constantes de integración:


[math]K=\frac{p1-p2}{\mu}[/math]
[math]C=0[/math]


Finalmente obtenemos la expresión f([math] \rho [/math]):


[math]f\left ( \rho \right )=\frac{1}{\mu }\left [ \frac{p2-p1}{4 }\rho ^{2} +p1-p2\right ][/math]


Para comprobar la condición de incomprensibilidad (el agua siempre ocupa el mismo volumen), vamos a comprobar que la divergencia del campo es nula.


[math]\bigtriangledown \cdot \overrightarrow{u}=0=\frac{1}{\rho } \left\{ \frac{\partial }{\partial \rho }\left ( \rho u_{\rho } \right )+ \frac{\partial }{\partial \theta }(u_{\theta })+ \frac{\partial }{\partial z}(\rho u_{z})\right\}[/math]


En conclusión, como f(p) no depende de "z", la divergencia es nula.

4 Campo de presiones y campo de velocidades

Vamos a dar como dato: [math] p1=4 [/math], [math] p2=1 [/math] y [math] \mu=1 [/math].


4.1 Campo de velocidades

Anteriormente hemos obtenido la función de velocidad. Vamos a proceder a sustituir los valores dados:´

[math]\vec{u}(\rho,\theta,z)=(\frac{-3\rho ^{2}}{4} +3)\vec{e_z}[/math]

Como podemos comprobar, a medida que nos acercamos al borde de la tubería, entra menos velocidad mientras que en el centro, hay mayor flujo de concentración de velocidades.

Campo de velocidades.
x=0:0.2:2; %Intervalo de 'ρ'
z=0:0.2:10; %Intervalo de 'z'
[X,Z]=meshgrid(x,z);
ux=(-3./4).*X.^2+3; %Dar valor a 'u'
uz=0.*zz;
hold on
quiver(X,Z,ux,uz)  %Representación
axis([0,5,0,12])
xlabel('ρ');
ylabel('z');
hold off
view(2)
title('Campo de velocidades')


4.2 Campo de presiones

La expresión del campo de presiones viene dada como:

[math] p(x,y,z)=p_1+\frac{(p_2-p_1)(z-1)}{2}[/math]

Sustituyendo los valores:

[math]p_1=4 \ {,} \ p_2=1 \ {,} \ μ=1[/math]

Obtenemos por tanto:

[math] p(x,y)=4-\frac{3z+3}{2}=\frac{11}{2}-\frac{3z}{2}[/math]


Como podemos comprobar, la presión no aumenta con el radio, sino con el parámetro "z".

Campo de presiones



z=0:0.1:10; %Altura 'z'
f=(-3)*z/2+(11/2); %Funcion del campo de presiones
plot(z,f,'r')
xlabel('Altura(z)');
ylabel('Presión(p)');
title(' Campo de presiones');


5 Líneas de corriente

Para dibujar las líneas de corriente del campo [math]\overrightarrow{u}[/math], debemos tener en cuenta que estas son tangenes a [math]\overrightarrow{u}[/math] en cada apunto.

Procedemos a calcular [math]\overrightarrow{v}[/math] que es ortogonal a [math]\overrightarrow{u}[/math] ya que se comprueba que: [math]\overrightarrow{v}=\overrightarrow{e_{\theta }}\cdot \overrightarrow{u}[/math].


Como hemos hallado anteriormente, la divergencia de [math]\overrightarrow{u}[/math] es nula. Consequentemente, el rotacional del campo [math]\overrightarrow{v}[/math] es nulo también.


La función de corriente o potencial escalar viene definido como: [math]\psi [/math],([math]\bigtriangledown \psi =\overrightarrow{v}[/math])

[math]\overrightarrow{v}=\begin{vmatrix} \overrightarrow{e_{\rho }} & \overrightarrow{e_{\theta }}&\overrightarrow{e_{z}} \\ 0&1 &0 \\ 0&0 &f\left ( \rho \right ) \end{vmatrix}=f\left ( \rho \right )\overrightarrow{e_{\rho }} \rightarrow \overrightarrow{v}=f\left ( \rho \right )\overrightarrow{e_{\rho }}[/math].


Sustituyendo [math]f\left ( \rho \right )[/math]:

[math]\overrightarrow{v}=\left ( \frac{p2-p1}{4\mu }\rho ^{2} +\frac{p1-p2}{\mu }\right )\overrightarrow{e_{\rho }}[/math].


[math]\bigtriangledown \psi=\frac{d\psi }{d\rho }\overrightarrow{e_{\rho }}+\frac{1}{\rho }\frac{d\psi }{d\theta }\overrightarrow{e_{\theta }}+\frac{d\psi }{dz}\overrightarrow{e_{z}}=\overrightarrow{v}[/math].


Obtenemos la relación:
[math]\frac{d\psi }{d\rho }=\frac{p2-p1}{4\mu }\rho ^{2} +\frac{p1-p2}{\mu }\ [/math].


Integrando el gradiente del potencial escalar:
[math]\psi=\frac{p2-p1}{12\mu }\rho ^{3} +\frac{p1-p2}{\mu }\rho [/math].


El resultado de la función potencial quedará:
[math] \psi= -\frac{1}{4} \rho ^{3} + 3\rho[/math].


Líneas de corriente
x=0:0.2:2; %Definimos 'ρ'
z=0:0.2:10; %Definimos 'z'
[X,Z]=meshgrid(x,z);
lineas=(-1./4).*X.^3+3.*X; %Definimos campo escalar  
contour(X,Z,lineas); 
axis([0,3,0,10]);
colorbar
title('Líneas de corriente');


6 Velocidad máxima del fluido y módulo de velocidad

Para obtener los máximos y mínimos de la velocidad vamos a optimizar la función de velocidad. Para ello, derivamos la misma y la igualamos a 0:


[math]\vec{u}(\rho,\theta,z)=(\frac{-3\rho ^{2}}{4 } +3)\vec{e_z}[/math]


donde:

[math] \frac{\partial \vec{u}} {\partial \rho}=\frac{-3\rho}{2}\vec{e_z}[/math]
[math] \frac{\partial \vec{u} }{\partial \theta}=0 [/math]
[math] \frac{\partial \vec{u} }{\partial z}=0 [/math]


Consequentemente, el punto donde la velocidad es máxima es en [math] \rho =0 [/math]

Como podemos observar en la gráfica la velocidad de nuestro fluido no es constante a lo largo de la tubería. De hecho, disminuye a medida que nos acercamos a los bordes de la tubería. Esto se debe a las fuerzas viscosas potenciadas por la diferencia de presión que debe mantener el fluido para ser viscoso.


Módulo velocidad de fluido



x=0:0.05:2; %Definimos rho
f=abs((-3/4)*x.^2+3); %Definimos funcion de velocidad
plot(x,f,'r');
title('Módulo de la velocidad')
xlabel('Radio')
ylabel('Velocidad')
axis([0,3,0,5])


7 Rotacional

Para hallar el rotacional, vamos a proceder con el cálculo:

[math]\nabla\times\vec u= \frac{1}{\rho}\begin{vmatrix} \vec{e_\rho} & \rho \vec{e_\theta} & \vec{e_z} \\ \frac{ \partial}{\partial \rho} & \frac{\partial }{\partial \theta} & \frac{\partial }{\partial z}\\ u_\rho & \rho u_\theta & u_z \end{vmatrix}[/math]

Sustituyendo en la fórmula:

[math]\nabla\times\vec u= \frac{1}{\rho}\begin{vmatrix} \vec{e_\rho} & \rho \vec{e_\theta} & \vec{e_z} \\ \frac{ \partial}{\partial \rho} & \frac{\partial }{\partial \theta} & \frac{\partial }{\partial z}\\ 0 & 0 & (\frac{-3\rho ^{2}}{4 } +3)\end{vmatrix}[/math]

Obtenemos que [math]\nabla\times\vec{u}=(\frac{3\rho}{2})\vec{e_\theta}[/math],


Para visualizar esto de manera gráfica, nuevamente a Octave,

Rotacional del campo velocidad
x=0:0.1:2; %Definimos 'rho'
 z=0:0.1:10; %Definimos 'z'
 [X,Z]=meshgrid(x,z);
 rot=abs(3.*X./2); %Calculamos el rotacional
 surf(X,Z,rot)
 colorbar
 view(2)
 axis([0,5,0,12])


8 Campo de temperaturas

En el enunciado nos viene dado la función Temperatura en coordenadas cilíndricas:

[math] T(ρ{,}θ{,}z)=log(1+ρ)e^{-(z-2)^2}[/math]


En la gráfica podemos observar como la temperatura es máxima para p=2, z=2; ya que los valores p y z son proporcionales a la temperatura.


Campo de temperaturas
p=0:0.01:2; %Definimos p
z=-2:0.05:10; %Definimos z
[P,Z]=meshgrid(p,z);
figure (1) 
a=log(1+P)*exp.^-(-Z-2).^2 %Funcion de temperatura
pcolor(Y,Z,p);
shading flat
grid on
axis([0,8,-1,10]);
10 colorbar 
11 figure(2)
12 contour(Y,Z,a,10,'k'); 
13 grid on
14 axis([0,8,-1,10]);


8.1 Curvas de nivel

Las curvas de nivel son líneas que concentran puntos con la misma temperatura. Se comprueba que tienen forma simétrica y logarítmica.


Curvas de nivel de campo de temperatura
x=0:0.05:2; %Definimos 'rho' 
z=0:0.05:1; %Definimos 'z'
[X,Z]=meshgrid(x,z);
figure (1)
p=log(1+P)*exp.^-(-Z-2).^2 %Funcion de temperatura
[pX,pZ]=gradient(p); %Funcion gradiente
hold on
quiver(X,Z,pX,pZ)
axis([0,8,-1,2]);
shading flat
grid on
hold off


9 Gradiente de temperatura

Teniendo la función de temperatura, calculamos su gradiente:

[math]\nabla T(ρ,θ,z)=(\frac{e^{-(z-2)^2}}{1+ρ)})\vec{e_\rho}-(2log(ρ+1)(z-2){e^{-(z-2)^2}})\vec{e_z}[/math]


La gráfica del gradiente de la temperatura, se ha representado utilizando el siguiente código:


derecha
rho=0:0.1:2; %Definimos 'rho'
z=0:0.1:10; %Definimos 'z'
[RHO,Z]=meshgrid(rho,z);
figure(1)
T=log(1+RHO).*exp(-(Z-2).^2); %Funcion temperatura
[TRHO,TZ]=gradient(T); %Funcion gradiente
hold on
quiver(RHO,Z,TRHO,TZ)
xlabel('rho');
ylabel('z');
axis([0,2,0,10]);
title('Gradiente de Temperatura')
shading flat
grid on
hold off


derecha
rho=0:0.1:2; %Definimos 'rho'
z=0:0.1:10; %Definimos 'z'
[RHO,Z]=meshgrid(rho,z);
figure(1)
T=log(1+RHO).*exp(-(Z-2).^2); %Definimos la funcion de temperatura
[TRHO,TZ]=gradient(T); %Definimos la funcion gradiente
hold on
quiver(RHO,Z,TRHO,TZ)
contour(RHO,Z,T,'k') %Colocamos las curvas de nivel superpuestas
xlabel('rho');
ylabel('z');
axis([0,2,0,10]);
title('Gradiente de Temperatura y Curvas de nivel')
shading flat
grid on
hold off



Como podemos verificar, las curvas de nivel son ortogonales al gradiente de la temperatura. Esto se debe a que el gradiente es el responsable de la dirección de crecimiento y decrecimiento del campo.

Por otro lado, podemos comprobar que si el gradiente es mayor, la temperatura tendrá una variación más rápida.

10 Caudal

El caudal es la cantidad de un fluido que atraviesa la tubería en un determinado tiempo. Para ello debemos tener en cuenta el flujo de volumen por unidad de tiempo. Vamos a representarlo de la siguiente forma:

[math]Q=\int_{S}\vec{v} \cdot d\vec{S}[/math]

donde [math]\vec{v}[/math] es el campo de velocidades.

El vector normal será perpendicular a la superficie por lo que hemos visto anteriormente de las curvas y gradientes.

El campo de velocidades finalmente será:

[math]\vec{u}(\rho,\theta,z)=(-\frac{3}{4}\rho ^{2} +3)\overrightarrow{e_{z}}[/math]
,

Finalmente resolviendo la integral doble obtendremos:

[math]Q\left ( m/s \right )=\int_{S}^{}\overrightarrow{u}\overrightarrow{n}dS=\int_{S}^{}\left (\frac{-3\rho^{2}}{4}+{3} \right )\overrightarrow{e_{z}}\cdot \overrightarrow{e_{z}} dS=\int_{0}^{2\Pi }\int_{0}^{2}\left ( \frac{-3\rho^{2}}{4}+{3} \right )d_{\rho }d_{\theta }=8{\Pi}=25.133\left ( m/s \right )[/math].