Fluido Alrededor de un Obstáculo Circular (Grupo 45)

Trabajo realizado por estudiantes
Título	Fluido Alrededor de un Obstáculo Circular (Grupo 45)
Asignatura	Teoría de Campos
Curso	2025-26
Autores	Jose Antonio Martín-Caro Xavier Grimalt Roig Uriel Hidalgo Marcos Emilio Tavío Pedro Comas Payeras
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

El objetivo de este artículo es estudiar el comportamiento de un fluido alrededor de un sólido circular.

Un fluido es una sustancia que puede deformarse continuamente bajo la aplicación de una fuerza de cizallamiento (es decir, una fuerza que actúa paralela a una superficie) sin mostrar resistencia permanente. A nivel físico, los fluidos pueden ser líquidos y gases, ya que ninguno de los dos puede conservar una forma estable. La diferencia entre ellos es que los primeros toman la forma del recipiente donde están, mientras que los segundos tienen tan poca unión entre sus partículas que pueden comprimirse y no tienen ni forma ni volumen propios.

1 .-Superficie Mallada

Se comienza realizando un mallado que describe los puntos interiores de la región ocupada por el fluido. Para llevar a cabo la representación de esta región se emplean coordenadas cilíndricas, definidas en el intervalo radial 1 ≤ r ≤ 5, que posteriormente se transforman a coordenadas cartesianas. Tras esta transformación, el dominio queda incluido en: (x,y) ∈ [−4,4] × [−4,4].

Mediante el siguiente código elaborado en MATLAB, se podrá visualizar la superficie de trabajo:

Figura 1 — Representación de Superficie Mallada

clear; clc;

%Mallado
r = linspace(1,5,60); %Radios entre 1 y 5
theta = linspace(0,2*pi,80); %Ángulos

% Mallado en coordenadas polares
[R,TH] = meshgrid(r,theta);

%Conversión a coordenadas cartesianas
X = R .* cos(TH);
Y = R .* sin(TH);
Z = zeros(size(R));
% Dibujar el mallado
figure; hold on;
mesh(X,Y,Z);

%Dibujar el círculo unidad
plot(1*cos(theta), 1*sin(theta), 'k', 'LineWidth', 2);

%Vista
axis equal;
axis([-4 4 -4 4]);
view(2); % Vista en planta
title('Mallado del Fluido (Región Exterior al Círculo Unidad)');
xlabel('X');
ylabel('Y');
hold off;

2 .-Función Potencial y Campo de Velocidades del Fluido.

Una vez definida la región donde se estudiará el comportamiento del fluido, se procede a analizar su dinámica a partir del campo de velocidades. Para ello, se utiliza una función potencial escalar: [math]\phi = \phi(\rho,\theta)[/math], definida en un dominio [math]D \subset \mathbb{R}^2[/math] expresado en coordenadas polares.

En este caso, la función potencial viene dada por: [math]\phi(\rho,\theta) = \left( \rho + \frac{1}{\rho} \right)\cos\theta[/math].

A partir de esta función se obtiene el campo de velocidades aplicando: [math]\vec{u} = \nabla\phi[/math], lo que permite describir la dirección y magnitud del movimiento del fluido, así como estudiar la relación entre las curvas de nivel del potencial y las líneas de corriente.

2.1 .-Representación de la Función Potencial

Para estudiar con mayor claridad la naturaleza del flujo, es útil examinar la forma que adopta la función potencial [math]\phi(\rho,\theta)[/math] dentro del dominio considerado. La representación gráfica de esta función permite identificar zonas donde el potencial crece o disminuye con mayor rapidez, así como patrones característicos que influyen en el comportamiento del fluido.

La construcción de estas gráficas se realiza mediante herramientas de visualización numérica, en este caso, MATLAB, que posibilitan generar superficies del potencial. Estas representaciones facilitan la interpretación del campo y sirven como apoyo previo al análisis del gradiente y de las velocidades resultantes.

Figura 2.1 — Curvas de nivel de la función potencial 𝜙 ( 𝜌 , 𝜃 )

clear; clc;

%Mallado en coordenadas polares
R = linspace(1,5,80); %Rho
T = linspace(0,2*pi,180); %Theta
[rho,theta] = meshgrid(R,T);

%Transformación a coordenadas cartesianas
X = rho.*cos(theta);
Y = rho.*sin(theta);

%Definición de la función potencial
phi = (rho+1./rho).*cos(theta);

%Representación de la función potencial (curvas de nivel)
figure;
contour(X, Y, phi, 60, 'LineWidth', 1.5);
hold on;
plot(1*cos(T), 1*sin(T), 'k', 'LineWidth', 2); %Círculo interior
axis equal;
colorbar;
title('Función Potencial');
xlabel('X');
ylabel('Y');
hold off;

2.2 .-Representación del Campo de Velocidades

Con la función potencial es posible determinar el campo de velocidades que describe el movimiento del fluido. Recordemos que la velocidad se calcula a partir del gradiente de la función potencial: [math]\vec{u} = \nabla\phi[/math].

Para ello se expresan las derivadas en coordenadas polares, donde el operador gradiente viene dado por: [math]\nabla\phi = \frac{\partial\phi}{\partial\rho}\,\vec{e}_\rho + \frac{1}{\rho}\frac{\partial\phi}{\partial\theta}\,\vec{e}_\theta[/math].

• La función potencial es: [math]\phi(\rho,\theta) = \left(\rho + \frac{1}{\rho}\right)\cos\theta[/math].

• Sus derivadas parciales son: [math]\frac{\partial\phi}{\partial\rho} = \left(1 - \frac{1}{\rho^{2}}\right)\cos\theta[/math]        [math]\frac{\partial\phi}{\partial\theta} = -\left(\rho + \frac{1}{\rho}\right)\sin\theta[/math].

Por tanto, las componentes del campo de velocidades en polares son: [math]u_{\rho} = \left(1 - \frac{1}{\rho^{2}}\right)\cos\theta[/math]         [math]u_{\theta} = -\left(1 + \frac{1}{\rho^{2}}\right)\sin\theta[/math].

Para representar el campo en MATLAB es más conveniente trabajar en coordenadas cartesianas. Usamos la transformación:

- [math]u_x = u_\rho\cos\theta - u_\theta\sin\theta[/math],

- [math]u_y = u_\rho\sin\theta + u_\theta\cos\theta[/math].

Tras realizar los cálculos:

- [math]u_x = (1 - \frac{1}{\rho^2})\cos^2\theta + (1 + \frac{1}{\rho^2})\sin^2\theta[/math],

- [math]u_y = -\frac{2}{\rho^2}\sin\theta\cos\theta[/math].

Figura 2.2 - Campo de velocidades alrededor del obstáculo circular

clear; clc;

%Mallado en polares
R = linspace(1,5,40);
T = linspace(0,2*pi,40);
[rho,theta] = meshgrid(R,T);

%Transformación a cartesianas
X = rho.*cos(theta);
Y = rho.*sin(theta);

%Función potencial
phi = (rho + 1./rho).*cos(theta);

%Componentes de la velocidad en polares
ur = (1 - 1./(rho.^2)).*cos(theta);
ut = -(1 + 1./(rho.^2)).*sin(theta);

%Conversión a componentes cartesianas
DX = ur.*cos(theta) - ut.*sin(theta);
DY = ur.*sin(theta) + ut.*cos(theta);

%Gráfica
figure;
contour(X, Y, phi, 40); hold on;
quiver(X, Y, DX, DY, 'k');
plot(cos(T), sin(T), 'r', 'LineWidth', 1.3); %Obstáculo circular
axis equal; axis([-5 5 -5 5]);
title('Campo de Velocidades');
xlabel('X'); ylabel('Y'); colorbar;
hold off;

Una característica interesante de este flujo es que las líneas de corriente coinciden con las trayectorias que seguirían partículas sin inercia, y estas son siempre perpendiculares a las curvas de nivel de la función potencial. Esto es una propiedad general de los flujos potenciales: el gradiente siempre apunta en la dirección de máxima variación del potencial, mientras que las curvas de nivel representan zonas donde el potencial es constante. Si se hace un zoom en cualquier área del diagrama se aprecia claramente que los vectores del campo de velocidades mantienen esta ortogonalidad en todo el dominio, especialmente alrededor del obstáculo circular donde el cambio de dirección es más brusco.

3 .-Rotacional y Divergencia

El rotacional y la divergencia de un campo vectorial proporcionan información esencial sobre el comportamiento físico del fluido que representan. La divergencia permite identificar si el fluido se comprime o se expande localmente, mientras que el rotacional muestra si las partículas experimentan algún tipo de giro o movimiento

3.1 .-Rotacional nulo

El rotacional muestra la dirección y la intensidad del giro del fluido en cada punto. Para analizar si el flujo induce rotación, se calcula el rotacional del campo de velocidades del siguiente modo:

• El campo de velocidades es: [math]\vec{u}=\nabla\phi=\frac{\partial\phi}{\partial\rho}\vec{e}_\rho+\frac{1}{\rho}\frac{\partial\phi}{\partial\theta}\vec{e}_\theta[/math] y sus componentes son: [math]u_\rho=\left(1-\frac{1}{\rho^2}\right)\cos\theta[/math], [math]u_\theta=-\left(1+\frac{1}{\rho^2}\right)\sin\theta[/math] y [math]u_z=0[/math].

• El rotacional en coordenadas cilíndricas es:[math]\nabla\times\vec{u}=\frac{1}{\rho}\begin{vmatrix}\vec e_\rho & \rho\vec e_\theta & \vec e_z \\[6pt] \frac{\partial}{\partial\rho} & \frac{\partial}{\partial\theta} & \frac{\partial}{\partial z} \\[6pt] \left(1-\frac{1}{\rho^2}\right)\cos\theta & -\rho\left(1+\frac{1}{\rho^2}\right)\sin\theta & 0 \end{vmatrix}=\frac{1}{\rho}\left[\vec e_\rho\cdot 0+\vec e_\theta\cdot 0+\vec e_z\left(\frac{\partial}{\partial\rho}\left(-\rho\left(1+\frac{1}{\rho^2}\right)\sin\theta\right)-\frac{\partial}{\partial\theta}\left(\left(1-\frac{1}{\rho^2}\right)\cos\theta\right)\right)\right]=0[/math]

Por lo tanto [math]\boxed { \nabla \times \bar { u } =0 } [/math] y esto significa que el campo es irrotacional por lo que las particulas del fluido no rotan sobre sí mismas al moverse por el flujo.

3.2 .-Divergencia nula

La divergencia de un campo vectorial es la cantidad que mide la diferencia entre el flujo que entra y el flujo que sale del volumen.

• [math]\nabla\cdot\vec{u}(\rho,\theta,z)=\frac{1}{\rho}\left\{\frac{\partial}{\partial\rho}(\rho u_\rho)+\frac{\partial}{\partial\theta}(u_\theta)+\frac{\partial}{\partial z}(\rho u_z)\right\}=\frac{1}{\rho}\left\{\frac{\partial}{\partial\rho}\left(\left(\rho-\frac{1}{\rho}\right)\cos\theta\right)+\frac{\partial}{\partial\theta}\left(-\left(1+\frac{1}{\rho^2}\right)\sin\theta\right)+0\right\}=\frac{1}{\rho}\left\{\left(1+\frac{1}{\rho^2}\right)\cos\theta-\left(1+\frac{1}{\rho^2}\right)\cos\theta\right\}=0[/math]

Como se observa, la divergencia resulta ser nula, es decir, [math] \boxed {\nabla \cdot \vec u = 0} [/math]. Esto implica que el volumen del fluido permanece constante: el flujo no se expande ni se contrae, por lo que el movimiento es incompresible.

4 .-Líneas de corriente del campo [math]\vec{u}[/math]

Las líneas de corriente de [math]\vec{u}[/math] son las curvas tangentes a los vectores velocidad de las partículas del fluido. Por lo que se comenzará calculando el campo [math]\vec{v}[/math], que en cada punto es ortogonal a [math]\boxed{\vec{u} : \vec{v} = \vec{k} \times \vec{u}}[/math] :

[math]\vec{v} =\begin{vmatrix} \vec e_\rho & \vec e_\theta & \vec e_z \\ 0 & 0 & 1 \\ (1-\frac{1}{\rho^2})\cos\theta & -(1+\frac{1}{\rho^2})\sin\theta & 0 \end{vmatrix} =(1+\tfrac{1}{\rho^2})\sin\theta\,\vec e_\rho +\left(1-\tfrac{1}{\rho^2}\right)\cos\theta\,\vec e_\theta [/math]

En el apartado 3.1 se ha comprobado que , el rotacional de [math]\vec{u}[/math] es nulo. Este resultado implica que [math]\vec{v}[/math] es irrotacional. Asimismo, [math]\vec{v}[/math] posee un potencial [math]\psi[/math], cuyo gradiente coincide con el propio campo.

Para determinar dicho potencial debe resolverse el sistema (se desprecia [math]V_z[/math]):

• [math]\frac{\partial \psi }{\partial \rho }=V_\rho [/math]

• [math]\frac{1}{\rho }\frac{\partial \psi}{\partial \theta}=V_\theta [/math]

Al haber despejado e integrado ambas ecuaciones [math]\frac{1}{\rho }[/math] se obtiene la función:[math] \boxed{ \psi=\sin\theta\left(\rho-\frac{1}{\rho}\right) } [/math]

Una vez calculado [math]\psi[/math] del campo [math]\vec{u}[/math], con la ayuda de MATLAB verificamos que los vectores de [math]\vec{v}[/math] son ortogonales a las líneas de corriente de [math]\psi[/math]. Se puede ver gráficamente, gracias a el siguiente codigo de matlab:

Figura 4 - Líneas de Corriente

clear;clc;
I=linspace(1,5,30);
J=linspace(0,2*pi,30);
[rho,theta]=meshgrid(I,J);
X=rho.*cos(theta);
Y=rho.*sin(theta);
%Función PSI
C=@(rho,theta)((rho-1./rho).*sin(theta));
Z=C(rho,theta);
contour(X,Y,Z,30);
hold on
%Gradiente de PSI (campo v)
DX= ((1+1./(rho.^2)).*sin(theta));
DY= ((1-1./(rho.^2)).*cos(theta));
quiver(X,Y,DX,DY);
plot(1.*cos(J),1.*sin(J),'k','lineWidth',2);
axis([-5,5,-5,5]);
colorbar;
title('Gráfica Líneas de Corriente');
xlabel('X');
ylabel('Y');
axis equal
view(2);
hold off

Aun cuando los vectores de [math]\vec{v}[/math] se presentan como ortogonales respecto a las líneas de corriente, los vectores de [math]\vec{u}[/math] deben ser tangentes a dichas curvas, manteniendo, a su vez, una perpendicularidad rigurosa con los vectores anteriormente citados de [math]\vec{v}[/math].Para poder demostrar esta afirmación graficamente, nos apoyamos en un codgio MATLAB.

Comparación entre [math]\vec v[/math] y [math]\vec u[/math]

clear;clc;
I=linspace(1,5,30);
J=linspace(0,2*pi,30);
[rho,theta]=meshgrid(I,J);
X=rho.*cos(theta);
Y=rho.*sin(theta);
%Función PSI
C=@(rho,theta)((rho-1./rho).*sin(theta));
Z=C(rho,theta);
contour(X,Y,Z,30);
hold on
%Gradiente de PSI (campo v)
DX= ((1+1./(rho.^2)).*sin(theta));
DY= ((1-1./(rho.^2)).*cos(theta));
quiver(X,Y,DX,DY);
%Gradiente de la función potencial phi (campo u)
DXX= (1-1./(rho.^2)).*cos(theta);
DYY= -(1+1./(rho.^2)).*sin(theta);
quiver(X,Y,DXX,DYY);
plot(1.*cos(J),1.*sin(J),'k','lineWidth',2);
view(2);
axis([-5,5,-5,5]);
colorbar;
title('Gráfica de la Función Potencial');
xlabel('X');
ylabel('Y');
axis equal;
hold off

Centro

Con el objetivo de una apreciación más detallada de la tangencialidad que forma [math]\vec{u}[/math] con las líneas de corriente y de ortogonalidad formada por [math]\vec{v}[/math] también con las las líneas de corriente, se ha ampliado la representación gráfica en un punto particular. Cabe destacar que esta relación se mantiene de manera uniforme en todo el campo, como se comprueba en la siguiente gráfica:

5 .-Puntos de Frontera S y Remanso

En este apartado se analizará la frontera [math]S[/math] del obstáculo circular de radio unidad.Se determinarán los puntos donde el módulo de la velocidad es máximo y mínimo.
Asimismo, se identificarán los puntos en los que la velocidad es nula, conocidos como puntos de remanso. Finalmente, se representarán gráficamente los puntos de remanso sobre el borde del obstáculo para comprobar su ubicación en el flujo.

5.1 .-Puntos de la frontera S

Sobre la frontera se cumple: [math]\rho = 1 . [/math]

Las componentes del campo de velocidades (ya calculadas previamente) eran: [math] u_\rho = \left( 1 - \frac{1}{\rho^2} \right)\cos\theta, \qquad u_\theta = -\left( 1 + \frac{1}{\rho^2} \right)\sin\theta. [/math]

En la frontera [math] \rho = 1 . [/math]

[math] u_{\rho}\big|_{\rho=1} = 0, \qquad u_{\theta}\big|_{\rho=1} = -2\sin\theta. [/math]

La rapidez sobre 𝑆 es:

[math] |\vec{u}| = 2\,|\sin\theta|. [/math]

• Máxima velocidad en 𝑆:

Ocurre cuando:

[math] |\sin\theta| = 1 \quad\Rightarrow\quad \theta = \frac{\pi}{2},\ \frac{3\pi}{2}. [/math]       En cartesianas: [math] (\,x,y\,)_{\text{max}} = (0,\,1),(0,\,-1). [/math]

• Mínima velocidad en 𝑆:

Ocurre cuando:

[math] |\sin\theta| = 0 \quad\Rightarrow\quad \theta = 0,\ \pi. [/math]           En cartesianas: [math] (\,x,y\,)_{\text{min}} = (1,\,0),(-1,\,0). [/math]

También se podría haber derivado la expresión de la rapidez e igualado a cero para obtener los extremos, pero ambos procedimientos —derivar o razonarlo directamente a partir de 2∣sin𝜃∣— conducen exactamente al mismo resultado.

5.2 .-Puntos de remanso

Los puntos de remanso sobre 𝑆 se obtienen imponiendo: [math] u_{\rho}\big|_{\rho=1}=0,\quad u_{\theta}\big|_{\rho=1}=0 \;\Longrightarrow\; -2\sin\theta=0 \;\Longrightarrow\; \sin\theta=0. [/math]
De aquí: [math] \theta = 0,\quad \theta=\pi. [/math]

En coordenadas cartesianas:

• [math]\theta = 0 \Rightarrow (x,y) = (1,0)[/math]
  
• [math]\theta = \pi \Rightarrow (x,y) = (-1,0)[/math]

Observaciones físicas:

En estos puntos la velocidad del fluido es nula respecto al sólido (puntos de estancamiento en la superficie del cilindro).
Según la ecuación de Bernoulli (en flujo potencial incompresible), un punto de remanso corresponde localmente a un máximo de presión estática.

Se pueden observar mejor estos puntos calculados visualmente a través de MATLAB:

Figura 5 - Visualización de los puntos de remanso y de velocidades máxima y mínima

clear; clc;

%Mallado en polares
R = linspace(1,5,40);
T = linspace(0,2*pi,40);
[rho,theta] = meshgrid(R,T);

%Transformación a cartesianas
X = rho.*cos(theta);
Y = rho.*sin(theta);

%Función potencial
phi = (rho + 1./rho).*cos(theta);

%Componentes de la velocidad en polares
ur = (1 - 1./(rho.^2)).*cos(theta);
ut = -(1 + 1./(rho.^2)).*sin(theta);

%Conversión a componentes cartesianas
DX = ur.*cos(theta) - ut.*sin(theta);
DY = ur.*sin(theta) + ut.*cos(theta);

%%PUNTOS DE LA FRONTERA S
%Velocidad máxima: (0,1) y (0,-1)
p_max = [0 1;
0 -1];

%Velocidad mínima: (1,0) y (-1,0)
p_min = [1 0;
-1 0];

%Gráfica
figure;
contour(X, Y, phi, 40); hold on;
quiver(X, Y, DX, DY, 'k');
plot(cos(T), sin(T), 'r', 'LineWidth', 1.3); %Obstáculo circular

%Dibujar puntos
plot(p_max(:,1), p_max(:,2), 'bo', 'MarkerFaceColor','b', 'MarkerSize',8);
plot(p_min(:,1), p_min(:,2), 'mo', 'MarkerFaceColor','m', 'MarkerSize',8);

legend('Líneas de potencial','Campo de velocidades',...
'Obstáculo','Vel. máxima','Vel. mínima');

axis equal; axis([-5 5 -5 5]);
title('Campo de Velocidades con Puntos Característicos');
xlabel('X'); ylabel('Y'); colorbar;
hold off;

6 .-Ecuación de Bernoulli

Con el fin de determinar los puntos donde el fluido alcanza mayor y menor presión, se considera una densidad constante igual a 2 (d=2). Además, debe cumplirse la ecuación de Bernoulli, de modo que [math]\tfrac{1}{2}d|\vec{u}|^2 + p = \text{cte}[/math] y para simplificar los cálculos, se asigna a dicha constante el valor 10.

Sustituyendo y despejando la presión, obtenemos finalmente: [math]p = 10 - |\vec{u}|^2[/math].

• El módulo del gradiente es: [math]|\vec{u}| = \sqrt{\left((1-\frac{1}{\rho^2})\cos\theta\right)^2 + \left(-(1+\frac{1}{\rho^2})\sin\theta\right)^2}[/math].

• El módulo del gradiente al cuadrado, es decir [math]|\vec{u}|^2[/math], resulta ser: [math]|\vec{u}|^2 = \cos^2\theta\left(1 - \frac{2}{\rho^2} + \frac{1}{\rho^4}\right) + \sin^2\theta\left(1 + \frac{2}{\rho^2} + \frac{1}{\rho^4}\right)[/math].

Por lo tanto, la presión queda dada por: [math]p = 10 - \left[\cos^2\theta\left(1 - \frac{2}{\rho^2} + \frac{1}{\rho^4}\right) + \sin^2\theta\left(1 + \frac{2}{\rho^2} + \frac{1}{\rho^4}\right)\right][/math].

Este campo de presiones se representa en la siguiente gráfica, obtenida mediante el código desarrollado en MATLAB.

Figura 6 - Campo de presiones

%% Parámetros del problema
Uinf = 1;            % Velocidad del flujo no perturbado
Rcil = 1;            % Radio del cilindro
NN = 420;            % Resolución de la malla

dens = 1.2;          % Densidad del fluido

%% Construcción de la malla 2D
xx = linspace(-3.2, 3.2, NN);
yy = linspace(-3.2, 3.2, NN);
[XX, YY] = meshgrid(xx, yy);

RR = hypot(XX, YY);
ANG = atan2(YY, XX);

%% Campo de velocidades
u_rad = Uinf .* cos(ANG) .* (1 - (Rcil^2 ./ RR.^2));
u_tan = -Uinf .* sin(ANG) .* (1 + (Rcil^2 ./ RR.^2));

Ux = u_rad .* cos(ANG) - u_tan .* sin(ANG);
Uy = u_rad .* sin(ANG) + u_tan .* cos(ANG);

Vnorm = sqrt(Ux.^2 + Uy.^2);
maskCil = (RR < Rcil);
Vnorm(maskCil) = NaN;

%% Coeficiente de presión
Cp_field = 1 - (Vnorm ./ Uinf).^2;

%% Presión dimensional
pres = 0.5 * dens * Uinf^2 .* Cp_field;
pres(maskCil) = NaN;

%% Gráfica
figure('Color', [1 1 1]);
contourf(XX, YY, Cp_field, 55, 'LineStyle', 'none');
colormap(turbo);
colorbar;

hold on
thetaPlot = linspace(0, 2*pi, 350);
plot(Rcil*cos(thetaPlot), Rcil*sin(thetaPlot), 'w', 'LineWidth', 2);

axis equal
xlabel('Coordenada X');
ylabel('Coordenada Y');
title('Distribución del coeficiente de presión alrededor de un cilindro');

7 .-Trayectoria de la Partícula

• 
• 
• 

8 .-Paradoja de D´Alembert

Como p [math]\vec{n}[/math] es la fuerza que ejerce el fluido en cada punto de la frontera, al sumar la proyección de todas estas fuerzas sobre la dirección [math]\vec{i}[/math] la resultante es nula. Por lo que, el fluido no realizara ninguna fuerza sobre el obstáculo en la dirección horizontal.

Para poder demostrarlo, hay que resolver la siguiente integral: [math] \int_{0}^{2\pi} p \cdot \vec{n} \cdot \vec{i} , d\theta [/math]

• [math]p = 5 - 4 \sin^2\theta [/math]

• [math] \vec{n} = - \vec{e}_\rho [/math]

• [math]\vec {i} = 1\cdot\vec {i}[/math]

Para hallar el valor de este vector en cilíndricas se usará la matriz de cambio de base previamente utilizada pero traspuesta: [math] \begin{pmatrix} cos(\theta) & sin(\theta) & 0 \\ -sin(\theta) & cos(\theta) & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} = cos(\theta)\vec {e}_\rho-sin(\theta)\vec {e}_\theta[/math]

- Finalmente se llega a [math]\vec {n}\cdot\vec {i}=-cos(\theta)[/math] pudiéndose así realizar la integral: [math]\int_0^{2\pi} p \cdot \vec{n} \cdot \vec{i} , d\theta = \int_0^{2\pi} (5 - 4 \sin^2\theta) \cdot (-\cos\theta) , d\theta = \int_0^{2\pi} (-5 \cos\theta + 4 \cos\theta \sin^2\theta) , d\theta [/math]

Resolviendo término a término:

• [math]\int_0^{2\pi} -5 \cos\theta , d\theta = [-5 \sin\theta]_0^{2\pi} = 0 [/math]

• [math] \int_0^{2\pi} 4 \cos\theta \sin^2\theta , d\theta = [\frac{4}{3} \sin^3\theta]_0^{2\pi} = 0[/math]

Por lo tanto, la integral total es nula: [math]\int_{0}^{2\pi} p \cdot \vec{n} \cdot \vec{i} , d\theta = 0 [/math]

9 .-Reanálisis de los apartados 2,3 y 4

9.1 .-Función Potencial y Campo de Velocidades

Una vez definida la región donde se estudiará el comportamiento del fluido, se procede a analizar su dinámica a partir del campo de velocidades. Para ello, se utiliza una función potencial escalar: [math]\phi = \phi(\rho,\theta)[/math], definida en un dominio [math]D \subset \mathbb{R}^2[/math] expresado en coordenadas polares.

En este caso, la función potencial viene dada por: [math]\phi(\rho,\theta) = \left( \rho + \frac{1}{\rho} \right)\cos\theta + \frac{\theta}{4\pi}[/math].

A partir de esta función se obtiene el campo de velocidades aplicando: [math]\vec{u} = \nabla\phi[/math], lo que permite describir la dirección y magnitud del movimiento del fluido, así como estudiar la relación entre las curvas de nivel del potencial y las líneas de corriente.

clear; clc;
I = linspace(1,5,30);
J = linspace(0,2*pi,30);
[rho,theta] = meshgrid(I,J);
hold on
X = rho.*cos(theta);
Y = rho.*sin(theta);
f = @(rho,theta)((rho + 1./rho).*cos(theta) + theta./(4*pi));
contour(X,Y,f(rho,theta),50)
plot(1.*cos(J),1.*sin(J),'k','lineWidth',2);
view(2);
axis([-5,5,-5,5]);
colorbar;
title('Función Potencial');
xlabel('X');
ylabel('Y');
axis equal;
hold off

Con la función potencial es posible determinar el campo de velocidades que describe el movimiento del fluido, al calcularse la velocidad a partir del gradiente de la función potencial: [math]\vec{u} = \nabla\phi[/math]

Para ello se expresan las derivadas en coordenadas polares, donde el operador gradiente viene dado por: [math]\boxed{\nabla\phi = \frac{\partial\phi}{\partial\rho}\,\vec{e}_{\rho} + \frac{1}{\rho}\frac{\partial\phi}{\partial\theta}\,\vec{e}_{\theta}}[/math]

La función potencial es: [math]\phi(\rho,\theta) = \left(\rho + \frac{1}{\rho}\right)\cos\theta + \frac{\theta}{4\pi}[/math]

Sus derivadas parciales son:

• [math]\frac{\partial\phi}{\partial\rho} = \left(1 - \frac{1}{\rho^{2}}\right)\cos\theta[/math]

• [math]\frac{\partial\phi}{\partial\theta} = -\left(\rho + \frac{1}{\rho}\right)\sin\theta + \frac{1}{4\pi}[/math].

Por tanto, las componentes del campo de velocidades en polares son:

• [math]u_{\rho} = \left(1 - \frac{1}{\rho^{2}}\right)\cos\theta[/math]

• [math]u_{\theta} = -\left(1 + \frac{1}{\rho^{2}}\right)\sin\theta + \frac{1}{4\pi\rho}[/math].

Para representar el campo en MATLAB es más conveniente trabajar en coordenadas cartesianas. Usamos la transformación:

- [math]u_x = u_\rho\cos\theta - u_\theta\sin\theta[/math],

- [math]u_y = u_\rho\sin\theta + u_\theta\cos\theta[/math].

Tras realizar los cálculos:

- [math]u_x = \left(1 - \frac{1}{\rho^2}\right)\cos^2\theta + \left(1 + \frac{1}{\rho^2}\right)\sin^2\theta - \frac{1}{4\pi\rho}\sin\theta[/math],

- [math]u_y = -\frac{2}{\rho^2}\sin\theta\cos\theta + \frac{1}{4\pi\rho}\cos\theta[/math].

clear; clc;
I = linspace(1,5,30);
J = linspace(0,2*pi,30);
[rho,theta] = meshgrid(I,J);
X = rho.*cos(theta);
Y = rho.*sin(theta);
f = @(rho,theta)((rho + 1./rho).*cos(theta) + theta./(4*pi));
contour(X,Y,f(rho,theta),50)
hold on
dphi_r = (1 - 1./(rho.^2)).*cos(theta);
dphi_t = -(1 + 1./(rho.^2)).*sin(theta) + 1./(4*pi*rho);
DX = dphi_r.*cos(theta) - dphi_t.*sin(theta);
DY = dphi_r.*sin(theta) + dphi_t.*cos(theta);
quiver(X,Y,DX,DY)
plot(1.*cos(J),1.*sin(J),'k','lineWidth',2);
axis([-5,5,-5,5]);
colorbar;
title('Gráfica de Velocidades');
xlabel('X');
ylabel('Y');
axis equal;
view(2);
hold off

9.2 .-Rotacional y Divergencia

- [math]\boxed{\text{Rotacional}}[/math]: El rotacional indica la dirección y la velocidad del giro del campo en cada punto.

En coordenadas cilíndricas se calcula como: [math]\nabla\times\vec{u}=\frac{1}{\rho}\begin{vmatrix}\vec{e}_\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]

Las componentes del campo obtenidas del potencial son:

• [math]u_\rho=\left(1-\frac{1}{\rho^{2}}\right)\cos\theta[/math]

• [math]u_\theta=-\left(1+\frac{1}{\rho^{2}}\right)\sin\theta+\frac{1}{4\pi\rho}[/math]

• [math]u_z=0[/math]

Sustituyendo: [math]\nabla\times\vec{u}=\frac{1}{\rho}\left[\vec{e}_\rho(0)+\vec{e}_\theta(0)+\vec{e}_z\left(\frac{\partial}{\partial\rho}(\rho u_\theta)-\frac{\partial u_\rho}{\partial\theta}\right)\right][/math]

Cálculo de los términos:

• [math]\frac{\partial}{\partial\rho}(\rho u_\theta)=-(1-\tfrac{1}{\rho^2})\sin\theta-\frac{1}{4\pi\rho^2}[/math]

• [math]\frac{\partial u_\rho}{\partial\theta}=-(1-\tfrac{1}{\rho^2})\sin\theta[/math]

Entonces: [math]\frac{\partial}{\partial\rho}(\rho u_\theta)-\frac{\partial u_\rho}{\partial\theta}=-\frac{1}{4\pi\rho^2}[/math] y por tanto: [math]\nabla\times\vec{u}=-\frac{1}{4\pi\rho^{3}}\vec{e}_z[/math]

Lejos del cilindro ([math]\rho\gg1[/math]) el campo es prácticamente irrotacional.

- [math]\boxed{\text{Divergencia}}[/math]: La divergencia mide el flujo neto entrante o saliente.

En coordenadas cilíndricas: [math]\nabla\cdot\vec{u}(\rho,\theta,z)=\frac{1}{\rho}\left\{\frac{\partial}{\partial\rho}(\rho u_\rho)+\frac{\partial u_\theta}{\partial\theta}+\frac{\partial u_z}{\partial z}\right\}[/math]

Sustituyendo:

• [math]\rho u_\rho=(\rho-\tfrac{1}{\rho})\cos\theta[/math]

• [math]\frac{\partial}{\partial\rho}(\rho u_\rho)=(1+\tfrac{1}{\rho^2})\cos\theta[/math]

• [math]\frac{\partial u_\theta}{\partial\theta}=-(1+\tfrac{1}{\rho^2})\cos\theta[/math]

Entonces: [math]\nabla\cdot\vec{u}=\frac{1}{\rho}\left[(1+\tfrac{1}{\rho^{2}})\cos\theta-(1+\tfrac{1}{\rho^{2}})\cos\theta\right]=0[/math] y por tanto, el campo es incompresible: [math]\nabla\cdot\vec{u}=0[/math].

9.3 .-Líneas de corriente del campo [math]\vec{u}[/math]

Las líneas de corriente de [math]\vec{u}[/math] son las curvas tangentes a los vectores de velocidad de las partículas del fluido. Para hallarlas simplemente hay que calcular el campo ortogonal a [math]\vec{u}[/math] resolviendo el siguiente producto vectorial:

[math]\vec{v}=\vec{k}\times\vec{u}=\begin{vmatrix}\vec{e}_\rho&\vec{e}_\theta&\vec{e}_z\\0&0&1\\\frac{\partial\varphi}{\partial\rho}&\frac{1}{\rho}\frac{\partial\varphi}{\partial\theta}&0\end{vmatrix}=-\left(\frac{1}{\rho}\frac{\partial\varphi}{\partial\theta}\right)\vec{e}_\rho+\left(\frac{\partial\varphi}{\partial\rho}\right)\vec{e}_\theta[/math]

La función potencial dada es: [math] \varphi(\rho,\theta)=\left(\rho+\frac{1}{\rho}\right)\cos\theta+\frac{\theta}{4\pi}. [/math]

Calculamos sus derivadas parciales:

• [math] \frac{\partial\varphi}{\partial\rho} = \left(1-\frac{1}{\rho^{2}}\right)\cos\theta [/math]
• [math] \frac{1}{\rho}\frac{\partial\varphi}{\partial\theta} = -(1+\frac{1}{\rho^{2}})\sin\theta + \frac{1}{4\pi\rho}. [/math]

Sustituyendo en la expresión de [math]\vec{v}[/math] obtenemos:

- [math]V_\rho = \left(1+\frac{1}{\rho^{2}}\right)\sin\theta - \frac{1}{4\pi\rho}[/math]

- [math]V_\theta = \left(1-\frac{1}{\rho^{2}}\right)\cos\theta[/math]

Como se ha comprobado en el apartado anterior, el rotacional de [math]\vec{v}[/math] es nulo. Esto quiere decir que [math]\vec{v}[/math] es irrotacional. Asimismo, [math]\vec{v}[/math] tiene un potencial llamado [math]\psi[/math], gradiente del propio [math]\vec{v}[/math], y que será descifrado resolviendo el siguiente sistema de ecuaciones diferenciales (se desprecia [math]V_z[/math] ya que se está trabajando en el plano): [math]\frac{\partial\psi}{\partial\rho}=V_\rho \qquad,\frac{1}{\rho}\frac{\partial\psi}{\partial\theta}=V_\theta[/math]

Tras despejar [math]\frac{1}{\rho}[/math] e integrar ambas ecuaciones obtenemos la función: [math]\boxed{\psi(\rho,\theta)=\left(\rho-\frac{1}{\rho}\right)\sin\theta-\frac{1}{4\pi}\ln(\rho)}[/math]

Una vez se obtiene tanto [math]\psi[/math] como [math]\vec{v}[/math], se introducen en MATLAB, así comprobando que efectivamente los vectores que componen [math]\vec{v}[/math] son ortogonales a las líneas de corriente de [math]\psi[/math]. Para ello se ha empleado el siguiente código:

clear; clc;
I = linspace(1,5,30);
J = linspace(0,2*pi,30);
[rho,theta] = meshgrid(I,J);
X = rho.*cos(theta);
Y = rho.*sin(theta);
C = @(rho,theta)( sin(theta).*(rho - 1./rho) + (theta/(4*pi)).*log(rho) );
Z = C(rho,theta);
contour(X,Y,Z,30)
hold on
dpsi_r = sin(theta).*(1 + 1./(rho.^2)) + (theta./(4*pi)).*(1./rho);
dpsi_t = cos(theta).*(rho - 1./rho) + (1/(4*pi))*log(rho);
DX = dpsi_r.*cos(theta) - (dpsi_t./rho).*sin(theta);
DY = dpsi_r.*sin(theta) + (dpsi_t./rho).*cos(theta);
quiver(X,Y,DX,DY)
plot(1.*cos(J),1.*sin(J),'k','lineWidth',2);
axis([-5,5,-5,5]);
colorbar;
title('Gráfica Líneas de Corriente');
xlabel('X');
ylabel('Y');
axis equal;
view(2);
hold off

A la vez que los vectores de [math]\vec{v}[/math] son ortogonales a las líneas de corriente, los de [math]\vec{u}[/math] deberían ser tangentes a estas, y a su vez perpendiculares a los ya mencionados vectores de [math]\vec{v}[/math]. Para demostrar esta afirmación gráficamente, se ha diseñado un nuevo código que permite observar los ángulos rectos que se forman:

clear; clc;
I = linspace(1,5,30);
J = linspace(0,2*pi,30);
[rho,theta] = meshgrid(I,J);
C = @(rho,theta)( sin(theta).*(rho - 1./rho) + (theta/(4*pi)).*log(rho) );
Z = C(rho,theta);
X = rho.*cos(theta);
Y = rho.*sin(theta);
contour(X,Y,Z,30)
hold on
dpsi_r = sin(theta).*(1 + 1./(rho.^2)) + (theta./(4*pi))./(rho);
dpsi_t = cos(theta).*(rho - 1./rho) + (1/(4*pi))*log(rho);
DX = dpsi_r.*cos(theta) - (dpsi_t./rho).*sin(theta);
DY = dpsi_r.*sin(theta) + (dpsi_t./rho).*cos(theta);
quiver(X,Y,DX,DY)
dphi_r = (1 - 1./(rho.^2)).*cos(theta);
dphi_t = -(1 + 1./(rho.^2)).*sin(theta) + 1./(4*pi*rho);
DXX = dphi_r.*cos(theta) - dphi_t.*sin(theta);
DYY = dphi_r.*sin(theta) + dphi_t.*cos(theta);
quiver(X,Y,DXX,DYY)
plot(1.*cos(J),1.*sin(J),'k','lineWidth',2);
title('Comparación entre v y u');
xlabel('X');
ylabel('Y');
axis([-5,5,-5,5]);
colorbar;
axis equal;
view(2);
hold off