Diferencia entre revisiones de «Hello World»

Revisión actual del 13:34 6 dic 2022

Flujo de Couette entre dos tubos concéntricos. Vamos a considerar el flujo de un fluido incompresible a través de dos cilindros concéntricos de manera que el exterior se mueve con velocidad angular constante en sentido antihorario mientras que el interior está fijo. Si suponemos que ambos cilindros tienen su eje en OX3 y pintamos la sección transversal (x3 = 0) el cilindro exterior queda proyectado sobre la la circunferencia ρ = 2 y el interior sobre la circunferencia ρ = 1. La velocidad angular cilindro exterior es ω > 0

Contenido

1 Dibujo de la sección trasversal
2 Cálculo de la velocidad de las partículas
3 Graficación del campo de velocidades
4 Líneas de corriente del flujo de Couette
5 Velocidad en función del radio
6 Calculo y graficación del rotacional de [math] \vec{u} [/math]
7 Temperatura

1 Dibujo de la sección trasversal

En primer lugar, con el fin de visualizar el flujo, cortamos los dos cilindros con el plano [math]x_3=0[/math], de forma que resulta la siguiente sección trasversal.

Sección trasversal de los dos cilindros

Para hallar esta figura, hemos ejecutado en MatLab el siguiente programa.

h=30;                      % definicion del intervalo
 u=linspace(1,2,h);         % pertenencia del parametro u [1,2]
 v=linspace(0,2*pi,2*h);    % pertenencia del parametro  v [0,2*pi]
 [U,V]=meshgrid(u,v);        % Matrices de coordenadas de U y V
 figure(1)
 X=U.*cos(V);                % parametrizacion
 Y=U.*sin(V);
 mesh(X,Y,0*X);               % Dibujo de la matriz
 axis([-3,3,-3,3])           % Selección de los ejes del dibujo
 view(2)                     % Elección de perspectiva

2 Cálculo de la velocidad de las partículas

2.1 Definición del campo de velocidades

Conocemos por la física del problema, que la velocidad de las partículas viene dada por [math]\vec{u}(ρ, θ) = f(ρ) \vec{e_θ} [/math], además de que su presión (p) es constante. Por otro lado, el campo de velocidades tiene que cumplir la ecuación de Navier-Stokes estacionaria:

[math](\vec{u} · ∇)\vec{u} + ∇p = µ∆\vec{u} [/math]

2.2 Cálculo del gradiente de la presión y de la parte convectiva

Al ser la presión constante, su gradiente es nulo: ∇p=0; y si despreciamos la parte convectiva [math](\vec{u} · ∇)\vec{u}=0 [/math] tendremos:

[math]µ∆\vec{u} =\vec{0}[/math]

Siendo µ el coeficiente de viscosidad del fluido y [math]∆\vec{u}[/math] el laplaciano vectorial del campo de velocidades.

2.3 Cálculo del laplaciano vectorial del capo de velocidades

El cálculo del laplaciano vectorial en coordenadas cartesianas es relativamente sencillo de calcular, ya que:

[math]∆\vec{u} = ∆(u_1\vec{i} + u_2\vec{j}+ u_3\vec{k}) = ∆u_1\vec{i} + ∆u_2\vec{j}+ ∆u_3\vec{k}[/math]

Pero como el campo de velocidades está dado en coordenadas cilindricas, resulta más sencillo calcularlo con la siguiente fórmula:

[math]\Delta\vec{u} = \nabla(\nabla \cdot \vec{u}) - \nabla \times (\nabla \times \vec{u})[/math]

2.3.1 Gradiente de la divergencia

En el primer sumando, aparece el gradiente de la divergencia, luego, en primer lugar, calculamos la divergencia. A modo de hipótesis, la divergencia ha de ser nula, ya que el fluido es incompresible.

[math] ∇ · \vec{u} = \frac{1}{ρ}[\frac{ \partial}{\partial ρ}(0) + \frac{ \partial}{\partial θ}(f(ρ)) + \frac{ \partial}{\partial z}(0) = 0 [/math]

Por lo tanto, se confirma la hiótesis, y tenemos la certeza de que se trata de un fluido incompresible. Además de esto:

[math] ∇(∇ · \vec{u})= ∇(0)= \vec{0} [/math]

, lo que concluye que el primer término en el calculo del laplaciano es nulo.

2.3.2 Rotacional del campo de velocidades

[math](\nabla\times\vec{u}) = \frac{1}{ρ}\begin{vmatrix} \vec{e_ρ} & ρ\vec{e_θ} & \vec{e_z} \\ \frac{ \partial}{\partial ρ} & \frac{\partial }{\partial θ} & \frac{\partial }{\partial z}\\ 0 & ρf(ρ) & 0 \end{vmatrix} = \frac{1}{ρ}\frac{ \partial (ρf(ρ)) }{\partial ρ}\vec{e_z} = [\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}]\vec{e_z}[/math]

2.3.3 Rotacional del rotacional del campo de velocidades

Una vez hallado el rotacional del campo de velocidades, seguimos el procedimiento calculando el rotacional del vector hallado:

[math]\nabla\times(\nabla\times\vec{u})= \frac{1}{ρ}\begin{vmatrix} \vec{e_ρ} & ρ\vec{e_θ} & \vec{e_z} \\ \frac{ \partial}{\partial ρ} & \frac{\partial }{\partial θ} & \frac{\partial }{\partial z}\\ 0 & 0 & [\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}] \end{vmatrix} = -\frac{1}{ρ}\frac{ \partial }{\partial ρ}(\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}) ρ\vec{e_θ} = \frac{1}{ρ}[\frac{f(ρ)}{ρ}-\frac{\partial(f(ρ))}{\partial ρ}-ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ}[/math]

2.3.4 Calculo final de Laplaciano vectorial

Una vez hallados todos los sumandos, sustituimos en la definición de laplaciano vectorial

[math]\Delta\vec{u} = \nabla(\nabla \cdot \vec{u}) - \nabla \times (\nabla \times \vec{u})[/math] [math]\Delta\vec{u} = \vec{0} - \frac{1}{ρ}[\frac{f(ρ)}{ρ}-\frac{\partial(f(ρ))}{\partial ρ}-ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ}[/math] [math]\Delta\vec{u} = \frac{1}{ρ}[-\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}+ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ}[/math]

2.4 Calculo final de la ecuación de Navier-Stocks

Conocidas ya todas los partes de la ecuación inicial, procedemos a resolverla, para hallar así el campo de velocidades:

[math](\vec{u} · ∇)\vec{u} + ∇p = µ∆\vec{u} [/math] [math] \vec + \vec{0} + {0}µ∆\vec{u} = \vec{0}[/math] [math]∆\vec{u} = \vec{0}[/math] [math]\frac{1}{ρ}[-\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}+ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ} = \vec{0} [/math] [math][-\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}+ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ} = \vec{0} [/math]

2.4.1 Obtención de una ecuación diferencial

Si tratamos de reducir la ecuación diferencial obtenida, llegamos a la siguiente conclusión:

[math] [-\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}+ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ} = [-\frac{f(ρ)}{ρ}+\frac{\partial}{\partial ρ}(ρ\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ} = \vec{0} [/math]

Este paso es sencillo de comprobar, ya que si hacemos la derivada del segundo sumando, inmediatamente resulta en la expresión anterior:

[math]\frac{\partial}{\partial ρ}(ρ\frac{\partial(f(ρ))}{\partial ρ}) = \frac{\partial(f(ρ))}{\partial ρ}+ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})[/math]

Luego, reordenamos para conseguir una ecuación diferencial de segundo orden

[math]\frac{f(ρ)}{ρ}-\frac{\partial}{\partial ρ}(ρ\frac{\partial(f(ρ))}{\partial ρ}) = 0 [/math] [math]\frac{f(ρ)}{ρ} = \frac{\partial}{\partial ρ}(ρ\frac{\partial(f(ρ))}{\partial ρ}) [/math]

2.4.2 Comprobación de una solución conocida

Tras despejar y reordenar los términos, resulta una compleja ecuación diferencial de la cual conocemos una solución:

[math]f(ρ) = aρ +\frac{b}{ρ},\hspace{20pt}a,b \in \mathbb{R} [/math]

Para comprobar que esta solución es valida, es necesario derivar la expresión de forma que dichas derivadas aparezcan en la ecuación diferencial obtenida.

[math]\frac{\partial f(ρ)}{\partial ρ} = a -\frac{b}{ρ^2},\hspace{20pt}a,b \in \mathbb{R} [/math] [math]\frac{\partial }{\partial ρ}(ρ\frac{\partial f(ρ)}{\partial ρ}) = \frac{\partial }{\partial ρ}(ρ(a -\frac{b}{ρ^2})) = \frac{\partial }{\partial ρ}(aρ -\frac{b}{ρ})) = a +\frac{b}{ρ^2},\hspace{20pt}a,b \in \mathbb{R}[/math]

Una vez halladas dichas derivadas, las introducimos en la ecuación y verificamos que:

[math] \frac{1}{ρ}(aρ +\frac{b}{ρ}) = a +\frac{b}{ρ^2} \Longrightarrow \frac{f(ρ)}{ρ} = \frac{\partial}{\partial ρ}(ρ\frac{\partial(f(ρ))}{\partial ρ}),\hspace{20pt}a,b \in \mathbb{R} [/math]

2.4.3 Búsqueda de las constantes "a" y "b"

Por ultimo, para hallar el campo de velocidades concreto, tenemos que imponer las condiciones que conocemos, ya que al resolver la ecuación diferencial de segundo orden han aparecido dos constantes desconocidas.

Condición 1: Si [math]ρ=1 \Longrightarrow \vec{u}=\vec{0} \Longrightarrow f(1)\vec{e_θ}=\vec{0} \Longrightarrow (a+b)\vec{e_θ}=\vec{0}\Longrightarrow a+b=0 [/math]

Condición 2: Si [math]ρ=2 \Longrightarrow \vec{u}= \omega \vec{e_θ} \Longrightarrow f(2)\vec{e_θ}=\omega \vec{e_θ} \Longrightarrow (2a+\frac{b}{2})\vec{e_θ}=\vec{0}\Longrightarrow 2a+\frac{b}{2}= 2\omega [/math]

Una vez impuestas estas dos condiciones, resulta un sistema de dos ecuaciones con dos incógnitas donde

[math] a = \frac{4}{3}\omega [/math] [math] b = -\frac{4}{3}\omega[/math]

concluyendo por lo tanto que:

[math]f(ρ) = \frac{4}{3}\omega ρ -\frac{4}{3ρ}[/math]

3 Graficación del campo de velocidades

Suponiendo que tanto la velocidad angular "ω" como el módulo de viscosidad "µ" son unitarios, dibujamos el campo de velocidades.

Campo vectorial de las velocidades del fluido

Para dibujar dicho gráfico, hemos introducido este código en el programa MatLab

w=1;                                      % Velocidad angular dada
 a=4/3*w;                                  % Definicion del parámetro a en función de la velocidad angular
 b=-a;                                     % Definicion del parámetro b en función de la velocidad angular
 n=30;                                     % Longitud de los intervalos
 u=linspace(1,2,n);                        % Definicioón del parametro u (rho)
 v=linspace(0,2*pi,n);                     % Definición del parametro v (theta)
 [U,V] = meshgrid(u,v);                    % Matriz de los parámetros
 x = U.*cos(V);                            % Coordenada x (cartesiana)
 y = U.*sin(V);                            % Coordenada y (cartesiana)
 fx = -sin(V)*((a)*U+b*(1./U));            % Campo vectorial en la dirección de x
 fy = cos(V)*((a)*U+b*(1./U));             % Campo vectorial en la dirección de y
 quiver(x,y,fx,fy);                        % Dibujo del campo vectorial

4 Líneas de corriente del flujo de Couette

Procedemos a calcular las líneas de corriente del campo [math]\vec{u} [/math], es decir las líneas que son tangentes a [math]\vec{u} [/math] en cada punto. Para ello, calculamos el vector [math]\vec{v} = \vec{e_z} \times \vec{u} [/math].

[math] \vec{v} = \vec{e_z} \times f(ρ)\vec{e_θ} =(\frac{4}{3}ρ- \frac{4}{3ρ})( \vec{e_z} \times \vec{e_θ}) = ( \frac{4}{3ρ}-\frac{4}{3}ρ)\vec{e_ρ} [/math]

4.1 irrotacionalidad de [math] \vec{v}[/math]

Calculamos el rotacional con el fin de confirmar que este es nulo:

[math](\nabla\times\vec{v}) = \frac{1}{ρ}\begin{vmatrix} \vec{e_ρ} & ρ\vec{e_θ} & \vec{e_z} \\ \frac{ \partial}{\partial ρ} & \frac{\partial }{\partial θ} & \frac{\partial }{\partial z}\\ \frac{4}{3ρ}-\frac{4}{3}ρ & 0 & 0 \end{vmatrix} = \frac{\partial }{\partial z}( \frac{4}{3ρ}-\frac{4}{3}ρ)ρ\vec{e_θ} - \frac{\partial }{\partial θ}( \frac{4}{3ρ}-\frac{4}{3}ρ)\vec{e_z} = \vec{0} [/math]

Lo que confirma la hipótesis: el rotacional de [math] \vec{v}[/math] es nulo.

4.2 Potencial escalar de [math] \vec{v}[/math]

Una vez hallado el vector [math]\vec{v} [/math], hallamos su potencial escalar, que representa la función de corriente de [math]\vec{u} [/math], es decir, hallamos la función ψ que cumpla [math]∇ψ = \frac{4}{3ρ}-\frac{4}{3}ρ\vec{v}[/math]

[math] ∇ψ = \vec{v} = (\frac{4}{3ρ} - \frac{4ρ}{3}) \vec{e_ρ} [/math] [math] ∇ψ = \frac{\partial ψ}{\partial ρ} \vec{e_ρ} + \frac{1}{ρ}\frac{\partial ψ}{\partial θ} \vec{e_θ} + \frac{\partial ψ}{\partial z}\vec{e_z} = (\frac{4}{3ρ} - \frac{4ρ}{3}) \vec{e_ρ}[/math] [math]\frac{\partial ψ}{\partial ρ} = \frac{4}{3ρ} - \frac{4ρ}{3}[/math] [math]ψ=\int (\frac{4}{3ρ}-\frac{4ρ}{3}) \partial ρ [/math] [math]ψ=\frac{4}{3}(ln(ρ)-\frac{ρ^2}{2}) \partial ρ [/math]

4.3 Graficación de las curvas de corriente

A continuación dibujamos la líneas de corriente.

Lineas de corriente del fluido

Para hallar esta figura, hemos ejecutado en MatLab el siguiente programa.

h=30;                        % definicion del intervalo
u=linspace(1,2,h);           % pertenencia del parametro u [1,2]
v=linspace(0,2*pi,2*h);      % pertenencia del parametro  v [0,2*pi]
[U,V]=meshgrid(u,v);         % Matrices de coordenadas de U y V
figure(1)
f = 4/3*(log(U)+(U.^2)./2);  % Campo escalar
X=U.*cos(V);                 % Parametrización
Y=U.*sin(V);
axis([-3,3,-3,3])            % Elección de los ejes
view(2)                      % Ver la figura desde arriba
contour(X,Y,f,10);           % Dibujo de las líneas de nivel
axis([-3,3,-3,3])            % Definición de los ejes del dibujo
colormap ocean               % Cambio del mapa de colores para una mejor visualización

5 Velocidad en función del radio

La velocidad del fluido es máxima cuando este está en contacto con el cilindro exterior, y es mínima cuando está en contacto con el cilindro interior. La variación de la velocidad en función del radio, viene dada por el módulo del campo de velocidades.

[math]|\vec{u}| = |(\frac{4}{3}\omega ρ -\frac{4}{3ρ})\vec{e_\theta}|= \frac{4}{3}\omega ρ -\frac{4}{3ρ}[/math]

Una vez hallado dicho campo escalar, procedemos a dibujarlo.

Campo escalar del módulo de la velocidad del fluido

Esto ha sido dibujado en MatLab con el siguiente programa. Comprobamos que el dibujo es independiente de ω desde esta perspectiva. En este dibujo, hemos definido ω= w=1.

w=1;                                      % Velocidad angular dada
 a=4/3*w;                                  % Definicion del parámetro a en función de la velocidad angular
 b=-a;                                     % Definicion del parámetro b en función de la velocidad angular
 n=30;                                     % Longitud de los intervalos
 u=linspace(1,2,n);                        % Definicioón del parametro u (rho)
 v=linspace(0,2*pi,n);                     % Definición del parametro v (theta)
 [U,V] = meshgrid(u,v);                    % Matriz de los parámetros
 x = U.*cos(V);                            % Coordenada x (cartesiana)
 y = U.*sin(V);                            % Coordenada y (cartesiana)
 f = a*U+b./U                              % Módulo de la velocidad es función del radio
 surf(x,y,f);                              % Dibujo del campo escalar
 axis([-3,3,-3,3])                         % Definición de los ejes
 view(2);                                  % Vista en 2D

6 Calculo y graficación del rotacional de [math] \vec{u} [/math]

El calculo del rotacional de [math] \vec{u} [/math] lo hemos calculado previamente para conseguir el laplaciano. Concluyendo que:

[math](\nabla\times\vec{u}) = [\frac{\frac{4}{3}\omega ρ -\frac{4}{3ρ}}{ρ}+\frac{\partial(\frac{4}{3}\omega ρ -\frac{4}{3ρ})}{\partial ρ}]\vec{e_z} = \frac{8}{3}\omega\vec{e_z}[/math]

Calculamos la norma del rotacional para saber cuales son los puntos con mayor rotacional.

[math] |\nabla\times\vec{u}| = |\frac{8}{3}\vec{e_z}| = \frac{8}{3}\omega\vec{e_z} [/math]

Como [math] |\nabla\times\vec{u}| = cte \hspace{5pt} \forall ρ,θ [/math]. Todos los puntos tienen igual rotacional. Por lo tanto, el campo del rotacional es constante y se ve de la siguiente forma:

Campo vectorial del rotacional de la velocidad

Este gráfico tiene una perspectiva diferente al resto ya que el rotacional es perpendicular al plano XOY, y por lo tanto, con la perspectiva anterior, el campo no se ve.

w=1;                                      % Velocidad angular dada
 n=30;                                     % Longitud de los intervalos
 u=linspace(1,2,n);                        % Definición del parámetro u (rho)
 v=linspace(0,2*pi,n);                     % Definición del parámetro v (theta)
 [U,V] = meshgrid(u,v);                    % Matriz de los parámetros
 x = U.*cos(V);                            % Coordenada x (cartesiana)
 y = U.*sin(V);                            % Coordenada y (cartesiana)
 z = zeros(n);
 fz = ones(n).*(8*w/3);                    % Campo vectorial en la dirección de y
 quiver3(x,y,z,zeros(n),zeros(n),fz);      % Dibujo del campo vectorial
 axis([-3,3,-3,3,-0.5,1])                  %Definición de los ejes del dibujo

7 Temperatura

Conocemos que la temperatura del sistema viene dada por la siguiente expresión [math]T(ρ,θ) = 1+ρ^2 sin^2θ e^−(ρ−3/2)2 [/math]

@@ Línea 6: / Línea 6: @@
 En primer lugar, con el fin de visualizar el flujo, cortamos los dos cilindros con el plano <math>x_3=0</math>, de forma que resulta la siguiente sección trasversal.
-[[Archivo:SeccionTrasversal.jpeg|sinmarco|centro]]
+[[Archivo:SeccionTrasversal.jpeg|400px|miniaturadeimagen|centro|Sección trasversal de los  dos cilindros]]
@@ Línea 64: / Línea 64: @@
 Una vez hallado el rotacional del campo de velocidades, seguimos el procedimiento calculando el rotacional del vector hallado:
-<center>  <math>\nabla\times(\nabla\times\vec{u})= \frac{1}{ρ}\begin{vmatrix} \vec{e_ρ} & ρ\vec{e_θ} & \vec{e_z} \\ \frac{ \partial}{\partial ρ} & \frac{\partial }{\partial θ} & \frac{\partial }{\partial z}\\ 0 & 0 & [\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}] \end{vmatrix} = -\frac{1}{ρ}ρ\frac{ \partial [\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}) }{\partial ρ}\vec{e_θ} = [\frac{f(ρ)}{ρ}-\frac{\partial(f(ρ))}{\partial ρ}-\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ}</math></center>
+<center>  <math>\nabla\times(\nabla\times\vec{u})= \frac{1}{ρ}\begin{vmatrix} \vec{e_ρ} & ρ\vec{e_θ} & \vec{e_z} \\ \frac{ \partial}{\partial ρ} & \frac{\partial }{\partial θ} & \frac{\partial }{\partial z}\\ 0 & 0 & [\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}] \end{vmatrix} = -\frac{1}{ρ}\frac{ \partial }{\partial ρ}(\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}) ρ\vec{e_θ} = \frac{1}{ρ}[\frac{f(ρ)}{ρ}-\frac{\partial(f(ρ))}{\partial ρ}-ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ}</math></center>
+====Calculo final de Laplaciano vectorial====
+Una vez hallados todos los sumandos, sustituimos en la definición de laplaciano vectorial
+<center><math>\Delta\vec{u} = \nabla(\nabla \cdot \vec{u}) - \nabla \times (\nabla \times \vec{u})</math></center>
+<center><math>\Delta\vec{u} = \vec{0} - \frac{1}{ρ}[\frac{f(ρ)}{ρ}-\frac{\partial(f(ρ))}{\partial ρ}-ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ}</math></center>
+<center><math>\Delta\vec{u} = \frac{1}{ρ}[-\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}+ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ}</math></center>
 ===Calculo final de la ecuación de Navier-Stocks===
-Conocidos ya todos los sumandos de la ecuación inicial, procedemos a resolverla:
+Conocidas ya todas los partes de la ecuación inicial, procedemos a resolverla, para hallar así el campo de velocidades:
 <center><math>(\vec{u} · ∇)\vec{u} + ∇p = µ∆\vec{u} </math></center>
-<center><math>µ∆\vec{u} = \vec{0}</math></center>
+<center><math> \vec + \vec{0} + {0}µ∆\vec{u} = \vec{0}</math></center>
 <center><math>∆\vec{u} = \vec{0}</math></center>
-<center><math>[\frac{f(ρ)}{ρ}-\frac{\partial(f(ρ))}{\partial ρ}-\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ} = \vec{0} </math></center>
+<center><math>\frac{1}{ρ}[-\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}+ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ} = \vec{0} </math></center>
-<center><math>\frac{f(ρ)}{ρ}-\frac{\partial}{\partial ρ}(ρ\frac{\partial(f(ρ))}{\partial ρ}) = \vec{0}</math></center>
+<center><math>[-\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}+ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ} = \vec{0} </math></center>
+====Obtención de una ecuación diferencial====
+Si tratamos de reducir la ecuación diferencial obtenida, llegamos a la siguiente conclusión:
+<center><math> [-\frac{f(ρ)}{ρ}+\frac{\partial(f(ρ))}{\partial ρ}+ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ} = [-\frac{f(ρ)}{ρ}+\frac{\partial}{\partial ρ}(ρ\frac{\partial(f(ρ))}{\partial ρ})]\vec{e_θ} = \vec{0} </math></center>
+Este paso es sencillo de comprobar, ya que si hacemos la derivada del segundo sumando, inmediatamente resulta en la expresión anterior:
+<center><math>\frac{\partial}{\partial ρ}(ρ\frac{\partial(f(ρ))}{\partial ρ}) = \frac{\partial(f(ρ))}{\partial ρ}+ρ\frac{\partial}{\partial ρ}(\frac{\partial(f(ρ))}{\partial ρ})</math></center>
+Luego, reordenamos para conseguir una ecuación diferencial de segundo orden
+<center><math>\frac{f(ρ)}{ρ}-\frac{\partial}{\partial ρ}(ρ\frac{\partial(f(ρ))}{\partial ρ}) = 0 </math></center>
 <center><math>\frac{f(ρ)}{ρ} = \frac{\partial}{\partial ρ}(ρ\frac{\partial(f(ρ))}{\partial ρ}) </math></center>
-===Cálculo final del campo de velocidades===
 ====Comprobación de una solución conocida====
 Tras despejar y reordenar los términos, resulta una compleja ecuación diferencial de la cual conocemos una solución:
@@ Línea 98: / Línea 114: @@
 a+b=0 </math>
-Condición 2: Si <math>ρ=2 \Longrightarrow \vec{u}= \omega \vec{e_θ} \Longrightarrow f(2)\vec{e_θ}=\omega \vec{e_θ} \Longrightarrow  (2a+\frac{b}{2})\vec{e_θ}=\vec{0}\Longrightarrow 2a+\frac{b}{2}=\omega </math>
+Condición 2: Si <math>ρ=2 \Longrightarrow \vec{u}= \omega \vec{e_θ} \Longrightarrow f(2)\vec{e_θ}=\omega \vec{e_θ} \Longrightarrow  (2a+\frac{b}{2})\vec{e_θ}=\vec{0}\Longrightarrow 2a+\frac{b}{2}= 2\omega </math>
 Una vez impuestas estas dos condiciones, resulta un sistema de dos ecuaciones con dos incógnitas donde
-<center><math>  a = \frac{2}{3}\omega </math></center>  <center><math> b = -\frac{2}{3}\omega</math></center>
+<center><math>  a = \frac{4}{3}\omega </math></center>  <center><math> b = -\frac{4}{3}\omega</math></center>
+concluyendo por lo tanto que:
+ <center><math>f(ρ) = \frac{4}{3}\omega ρ -\frac{4}{3ρ}</math></center>
 ==Graficación del campo de velocidades==
@@ Línea 108: / Línea 127: @@
 Suponiendo que tanto la velocidad angular "ω" como el módulo de viscosidad "µ" son unitarios, dibujamos el campo de velocidades.
-[[Archivo:CampoDeVelocidades2.jpg|sinmarco|centro]]
+[[Archivo:CampoDeVelocidades3.jpg|400px|miniaturadeimagen|centro|Campo vectorial de las velocidades del fluido]]
 Para dibujar dicho gráfico, hemos introducido este código en el programa MatLab
 {{matlab|codigo=
   w=1;                                      % Velocidad angular dada
-  a=2/3*w;                                  % Definicion del parámetro a en función de la velocidad angular
+  a=4/3*w;                                  % Definicion del parámetro a en función de la velocidad angular
   b=-a;                                     % Definicion del parámetro b en función de la velocidad angular
   n=30;                                     % Longitud de los intervalos
@@ Línea 126: / Línea 145: @@
 }}
-==Graficación de las lineas de corriente==
+==Líneas de corriente del flujo de Couette ==
-Procedemos a dibujar las líneas de corriente del campo <math>\vec{u} </math>, es decir las líneas que son tangentes
+Procedemos a calcular las líneas de corriente del campo <math>\vec{u} </math>, es decir las líneas que son tangentes
 a <math>\vec{u} </math> en cada punto. Para ello, calculamos el vector <math>\vec{v} = \vec{e_z} \times \vec{u} </math>.
-<center><math> \vec{v} = \vec{e_z} \times f(ρ)\vec{e_θ} =(\frac{2}{3}ρ- \frac{2}{3ρ})( \vec{e_z} \times \vec{e_θ}) = ( \frac{2}{3ρ}-\frac{2}{3}ρ)\vec{e_ρ}  </math></center>
+<center><math> \vec{v} = \vec{e_z} \times f(ρ)\vec{e_θ} =(\frac{4}{3}ρ- \frac{4}{3ρ})( \vec{e_z} \times \vec{e_θ}) = ( \frac{4}{3ρ}-\frac{4}{3}ρ)\vec{e_ρ}  </math></center>
-Una vez hallado el vector <math>\vec{v} </math>, hallamos su potencial escalar, que representa la función de corriente de <math>\vec{u} </math>, es decir, hallamos la función ψ que cumpla <math>∇ψ = \vec{v}</math>
+===irrotacionalidad de <math> \vec{v}</math>===
-<center><math> ∇ψ = \vec{v} = (\frac{2}{3ρ} - \frac{2ρ}{3}) \vec{e_ρ} </math></center>
+Calculamos el rotacional con el fin de confirmar que este es nulo:
+<center>  <math>(\nabla\times\vec{v}) = \frac{1}{ρ}\begin{vmatrix} \vec{e_ρ} & ρ\vec{e_θ} & \vec{e_z} \\ \frac{ \partial}{\partial ρ} & \frac{\partial }{\partial θ} & \frac{\partial }{\partial z}\\ \frac{4}{3ρ}-\frac{4}{3}ρ & 0 & 0 \end{vmatrix} = \frac{\partial }{\partial z}( \frac{4}{3ρ}-\frac{4}{3}ρ)ρ\vec{e_θ}  - \frac{\partial }{\partial θ}( \frac{4}{3ρ}-\frac{4}{3}ρ)\vec{e_z} = \vec{0} </math></center>
+Lo que confirma la hipótesis: el rotacional de <math> \vec{v}</math> es nulo.
+===Potencial escalar de <math> \vec{v}</math>===
+Una vez hallado el vector <math>\vec{v} </math>, hallamos su potencial escalar, que representa la función de corriente de <math>\vec{u} </math>, es decir, hallamos la función ψ que cumpla <math>∇ψ = \frac{4}{3ρ}-\frac{4}{3}ρ\vec{v}</math>
+<center><math> ∇ψ = \vec{v} = (\frac{4}{3ρ} - \frac{4ρ}{3}) \vec{e_ρ} </math></center>
+<center><math> ∇ψ = \frac{\partial ψ}{\partial ρ} \vec{e_ρ} + \frac{1}{ρ}\frac{\partial ψ}{\partial θ} \vec{e_θ} + \frac{\partial ψ}{\partial z}\vec{e_z}  =  (\frac{4}{3ρ} - \frac{4ρ}{3}) \vec{e_ρ}</math></center>
+<center><math>\frac{\partial ψ}{\partial ρ} = \frac{4}{3ρ} - \frac{4ρ}{3}</math></center>
+<center><math>ψ=\int (\frac{4}{3ρ}-\frac{4ρ}{3}) \partial ρ </math></center>
+<center><math>ψ=\frac{4}{3}(ln(ρ)-\frac{ρ^2}{2}) \partial ρ </math></center>
+=== Graficación de las curvas de corriente===
+A continuación dibujamos la líneas de corriente.
+[[Archivo:LineasDeCorriente2.jpg|400px|miniaturadeimagen|centro|Lineas de corriente del fluido]]
+Para hallar esta figura, hemos ejecutado en MatLab el siguiente programa.
+{{matlab|codigo=
+h=30;                        % definicion del intervalo
+u=linspace(1,2,h);           % pertenencia del parametro u [1,2]
+v=linspace(0,2*pi,2*h);      % pertenencia del parametro  v [0,2*pi]
+[U,V]=meshgrid(u,v);         % Matrices de coordenadas de U y V
+figure(1)
+f = 4/3*(log(U)+(U.^2)./2);  % Campo escalar
+X=U.*cos(V);                 % Parametrización
+Y=U.*sin(V);
+axis([-3,3,-3,3])            % Elección de los ejes
+view(2)                      % Ver la figura desde arriba
+contour(X,Y,f,10);           % Dibujo de las líneas de nivel
+axis([-3,3,-3,3])            % Definición de los ejes del dibujo
+colormap ocean               % Cambio del mapa de colores para una mejor visualización
+}}
+== Velocidad en función del radio==
+La velocidad del fluido es máxima cuando este está en contacto con el cilindro exterior, y es mínima cuando está en contacto con el cilindro interior. La variación de la velocidad en función del radio, viene dada por el módulo del campo de velocidades.
+ <center><math>|\vec{u}| = |(\frac{4}{3}\omega ρ -\frac{4}{3ρ})\vec{e_\theta}|= \frac{4}{3}\omega ρ -\frac{4}{3ρ}</math></center>
+Una vez hallado dicho campo escalar, procedemos a dibujarlo.
+[[Archivo:GRADVEL3.jpg|400px|miniaturadeimagen|centro|Campo escalar del módulo de la velocidad del fluido]]
+Esto ha sido dibujado en MatLab con el siguiente programa. Comprobamos que el dibujo es independiente de ω desde esta perspectiva. En este dibujo, hemos  definido ω= w=1.
+{{matlab|codigo=
+ w=1;                                      % Velocidad angular dada
+ a=4/3*w;                                  % Definicion del parámetro a en función de la velocidad angular
+ b=-a;                                     % Definicion del parámetro b en función de la velocidad angular
+ n=30;                                     % Longitud de los intervalos
+ u=linspace(1,2,n);                        % Definicioón del parametro u (rho)
+ v=linspace(0,2*pi,n);                     % Definición del parametro v (theta)
+ [U,V] = meshgrid(u,v);                    % Matriz de los parámetros
+ x = U.*cos(V);                            % Coordenada x (cartesiana)
+ y = U.*sin(V);                            % Coordenada y (cartesiana)
+ f = a*U+b./U                              % Módulo de la velocidad es función del radio
+ surf(x,y,f);                              % Dibujo del campo escalar
+ axis([-3,3,-3,3])                         % Definición de los ejes
+ view(2);                                  % Vista en 2D
+}}
+== Calculo y graficación del rotacional de <math> \vec{u} </math>==
+El calculo del rotacional de <math> \vec{u} </math> lo hemos calculado previamente para conseguir el laplaciano. Concluyendo que:
+<center>  <math>(\nabla\times\vec{u}) = [\frac{\frac{4}{3}\omega ρ -\frac{4}{3ρ}}{ρ}+\frac{\partial(\frac{4}{3}\omega ρ -\frac{4}{3ρ})}{\partial ρ}]\vec{e_z} = \frac{8}{3}\omega\vec{e_z}</math></center>
+Calculamos la norma del rotacional para saber cuales son los puntos con mayor rotacional.
+<center>  <math> |\nabla\times\vec{u}| = |\frac{8}{3}\vec{e_z}| = \frac{8}{3}\omega\vec{e_z} </math></center>
+Como <math> |\nabla\times\vec{u}| = cte  \hspace{5pt} \forall ρ,θ  </math>. Todos los puntos tienen igual rotacional. Por lo tanto, el campo del rotacional es constante y se ve de la siguiente forma:
+[[Archivo:CampoRotacional3.jpg|400px|miniaturadeimagen|centro|Campo vectorial del rotacional de la velocidad]]
+Este gráfico tiene una perspectiva diferente al resto ya que el rotacional es perpendicular al plano XOY, y por lo tanto, con la perspectiva anterior, el campo no se ve.
+{{matlab|codigo=
+ w=1;                                      % Velocidad angular dada
+ n=30;                                     % Longitud de los intervalos
+ u=linspace(1,2,n);                        % Definición del parámetro u (rho)
+ v=linspace(0,2*pi,n);                     % Definición del parámetro v (theta)
+ [U,V] = meshgrid(u,v);                    % Matriz de los parámetros
+ x = U.*cos(V);                            % Coordenada x (cartesiana)
+ y = U.*sin(V);                            % Coordenada y (cartesiana)
+ z = zeros(n);
+ fz = ones(n).*(8*w/3);                    % Campo vectorial en la dirección de y
+ quiver3(x,y,z,zeros(n),zeros(n),fz);      % Dibujo del campo vectorial
+ axis([-3,3,-3,3,-0.5,1])                  %Definición de los ejes del dibujo
+}}
-<center><math> ∇ψ = \frac{\partial ψ}{\partial ρ} \vec{e_ρ} + \frac{1}{ρ}\frac{\partial ψ}{\partial θ} \vec{e_θ} + \frac{\partial ψ}{\partial z}\vec{e_z} 0  =  (\frac{2}{3ρ} - \frac{2ρ}{3}) \vec{e_ρ}</math></center>
+== Temperatura ==
-<center><math>\frac{\partial ψ}{\partial ρ} = \frac{2}{3ρ} - \frac{2ρ}{3}</math></center>
+Conocemos que la temperatura del sistema viene dada por la siguiente expresión <math>T(ρ,θ) = 1+ρ^2 sin^2θ e^−(ρ−3/2)2 </math>
-<center><math>ψ=\int (\frac{2}{3ρ}-\frac{2ρ}{3}) \partial ρ </math></center>
-<center><math>ψ=\frac{2}{3}(ln(ρ)-\frac{ρ^2}{2}) \partial ρ </math></center>