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

Trabajo realizado por estudiantes
Título	Fluido alrededor de un Obstáculo Circular. Grupo 33-A
Asignatura	Teoría de Campos
Curso	2023-24
Autores	Beatrice Laval González Mario Quero González Maximiliano Rodríguez Ruiz Gerardo Rodríguez Socas
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

En este artículo se estudiará el comportamiento de un fluido alrededor de un sólido circular.
Se define un fluido como aquel medio continuo cuya propiedad definitoria es la posibilidad de cambiar de forma sin que existan fuerzas restitutivas que tiendan a recuperar la forma original. De este modo, el término fluido hace referencia a aquellos gases y líquidos que carecen de rigidez y elasticidad. Además, son idealmente incompresibles ya que su volumen no disminuye al aplicarle fuerzas. Esta última cualidad, la tomaremos como objeto de estudio.

1 .-Superficie Mallada

Se comenzará realizando un mallado que describa los puntos interiores de la región ocupada por el fluido. Siendo este externo al obstáculo interior con una longitud de radio unitario. Para la representación del mismo, se emplearán coordenadas cilíndricas definidas en el intervalo Estando a su vez definidas en cartesianas en el intervalo (x,y) [5,5] x [-5,5].
Mediante el siguiente código elaborado en Matlab, se podrá visualizar la superficie de trabajo:

Representación del Mallado

clear;clc;
%Se parametrizan los valores de las circunferencias
I=linspace(1,5,50);   
J=linspace(0,2*pi,50);
%Se define el mallado del fluido
[rho,theta]=meshgrid(I,J);  
hold on
%Se definen X, Y y Z en coordenadas cilíndricas y se genera el mallado
X=rho.*cos(theta);      
Y=rho.*sin(theta);
Z=0.*rho;
mesh(X,Y,Z);
%Se crea el obstáculo
plot(1*cos(theta),1*sin(theta),'k','lineWidth',1);  
axis([-5,5,-5,5]);  
view(2); 
title('Mallado del Fluido');
xlabel 'X'
ylabel 'Y'
axis equal
hold off

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

Tras haber definido la región de estudio, se analizará el comportamiento del fluido. Para ello, será necesario calcular el campo de velocidades. El dominio del campo escalar estará definido por la función potencial siendo este [math] D\subset \mathbb{R}^{3} [/math] . La función potencial con la que se trabajará para conocer su comportamiento será: [math] \varphi (\rho ,\theta)=(\rho +\frac{1}{\rho})\cos (\theta )-\theta [/math]

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

En primer lugar, se mostrará la representación de la función potencial. Todo esto será desarrollado mediante funciones empleadas con la ayuda del programa Matlab:

Lineas de Nivel de la Función Potencial

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);
%Se introduce la función potencial
f=@(x,y)(x+(1./x)).*cos(y)-y;
%Se generan las líneas de nivel de la función
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

2.2 .-Representación del Campo de Velocidades

Sabiendo que la velocidad de las partículas de este fluido es el gradiente de la función potencial tal que: [math] \nabla \varphi=\vec u= \frac{\partial \vec{u}}{\partial \rho }\vec e_\rho + \frac{1}{\rho }\frac{\partial \vec{u}}{\partial \theta }\vec e_\theta+\frac{\partial \vec{u}}{\partial z}\vec e_z = \left ( 1-\frac{1}{\rho^2} \right )\cdot\cos(\theta)\vec e_\rho - [\left ( 1+\frac{1}{\rho^2} \right )\cdot\sin(\theta)-\frac{1}{\rho}]\vec e_\theta [/math]

A continuación, se mostrará una representación de los vectores del campo de velocidades siendo estos ortogonales a las curvas de nivel de la función potencial. Todo esto será desarrollado mediante funciones empleadas con la ayuda del programa Matlab, para trabajar con este programa se expresa el gradiente en coordenadas cartesianas, consiguiendo este cambio con ayuda de la siguiente matriz de cambio de base: [math] \begin{pmatrix} cos(\theta) & - sin(\theta) & 0 \\ sin(\theta) & cos(\theta) & 0 \\ 0 & 0 & 1 \end{pmatrix}[/math]

Velocidad de las Partículas

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=@(x,y)((x+(1./x)).*cos(y))-y;
contour(X,Y,f(rho,theta),50);
hold on
%Gradiente de la función potencial
DX=((cos(theta).^2).*(1-1./(rho.^2)))+((sin(theta).^2).*(1+ 1./(rho.^2)))+(sin(theta)./rho);
DY=((sin(theta).*cos(theta)).*(1-1./(rho.^2)))+((sin(theta).*cos(theta)).*(-1-1./(rho.^2)))-(cos(theta)./rho);
quiver(X,Y,DX,DY);
plot(1*cos(J),1*sin(J),'k','lineWidth',1);
axis([-5,5,-5,5]);
colorbar;
title('Gráfica de la Función Potencial');
xlabel('X');
ylabel('Y');
axis equal
view(2);
hold off

Si se hiciera zoom en cualquier punto de la gráfica del campo de velocidades, se podría observar que los vectores que representan la velocidad de las partículas del fluido, son siempre ortogonales a las curvas de nivel de la función potencial a lo largo de todo el movimiento alrededor del obstáculo circular. Como se puede observar en la siguiente imagen:

centro

3 .-Vector Normal a la Frontera del Obstáculo Circular

Es posible verificar que el campo de velocidades [math] \nabla \varphi=\vec {u} [/math] es ortogonal a las curvas de nivel del campo escalar mediante el vector normal [math] \vec n [/math]. Siendo [math] \vec n [/math] un vector normal cualquiera a los puntos del obstáculo .Si [math] \vec u \cdot \vec n = 0 [/math] , implica que son ortogonales entre sí. Al hacer los cálculos en coordenadas cilíndricas, el vector ortogonal a las curvas de nivel será: -[math] \vec e_\rho [/math]

Suponiendo el gradiente en ρ=1 y multiplicándolo por el vector ortogonal anteriormente nombrado, se puede comprobar mediante un breve cálculo la condición de ortogonalidad:
[math](\vec u(\rho,\theta,z)=\left ( 1-\frac{1}{\rho^2} \right )\cdot\cos(\theta)\vec e_\rho - [\left ( 1+\frac{1}{\rho^2} \right )\cdot\sin(\theta)-\frac{1}{\rho}]\vec e_\theta) \cdot (-\vec e_\rho) = 0 [/math]

Si calculamos el valor de [math] \vec u [/math] en la frontera del obstáculo se obtendría: [math]\vec u(1,\theta,z)=\nabla \varphi(1,\theta,z)= (- 2\cdot\sin(\theta) - 1 )\vec e_\theta [/math]

4 .-Valor aproximado de [math]\vec{u}[/math] lejos del Obstáculo

Si se está lejos del obstáculo, [math] \rho [/math] aumenta de tal forma, que [math] \frac{1}{\rho} [/math] pasa a considerarse despreciable. Por tanto, la función potencial pasaría a ser la siguiente: [math] \varphi_2 (\rho ,\theta)=\rho\cos (\theta )-\theta [/math].
Al cambiar la función potencial cambia el gradiente: [math] \varphi_2= \vec u_2 =\frac{\partial \vec{u}}{\partial \rho }\vec e_\rho + \frac{1}{\rho }\frac{\partial \vec{u}}{\partial \theta }\vec e_\theta+\frac{\partial \vec{u}}{\partial z}\vec e_z= (cos(\theta))\vec e_\rho - (sin(\theta) + \frac{1}{\rho})\vec e_\theta [/math].

Para conocer el valor del gradiente en cartesianas, hay que hacer un cambio de base, de nuevo con ayuda de la siguiente matriz: [math] \begin{pmatrix} cos(\theta) & - sin(\theta) & 0 \\ sin(\theta) & cos(\theta) & 0 \\ 0 & 0 & 1 \end{pmatrix}[/math] ,

obteniendo los siguientes resultados: [math] \vec u (\vec i , \vec j, \vec k )= (cos^2(\theta) + sin^2(\theta) + \frac{sin(\theta)}{\rho}) \vec i + (sin(\theta)cos(\theta) - sin(\theta)cos(\theta) - \frac{cos(\theta)}{\rho}) \vec j + 0 \vec k = (1 + \frac{sin(\theta)}{\rho}) \vec i - (\frac{cos(\theta)}{\rho}) \vec j [/math]

5 .-Comprobación de Rotacional y Divergencia Nula

El rotacional y la divergencia de un campo vectorial proporcionan información sobre las características físicas del fluido al que estudia y define.

5.1 .-Rotacional

Por su parte, el rotacional indica la dirección y la velocidad del giro del campo en cada punto. En él se considera el movimiento del fluido y su tendencia a inducir rotación en cada punto. Se calcula el rotacional, tal que:

[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}=\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}\\ cos\theta (1-\frac{1}{\rho^2}) & -(1+\frac{1}{\rho^2}) \rho sin\theta + 1 & 0 \end{vmatrix} = \frac{1}{\rho }[ (\vec e_\rho \cdot 0) + (\vec e_\theta \cdot 0) + \vec e_z \cdot (-sin\theta(1-\frac{1}{\rho^2}) + sin\theta(1-\frac{1}{\rho^2}))] = 0 [/math]

Por tanto, se puede concluir que es un campo irrotacional , dado que [math]\boxed { \nabla \times \bar { u } =0 } [/math] ,y, por ende, que las partículas del fluido no giran.

5.2 .-Divergencia

La divergencia de un campo vectorial es una magnitud escalar que compara el flujo entrante y el flujo saliente en una superficie.
[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( (\rho - \frac{1}{\rho^2})\cos(\theta) \right) + \frac{\partial}{\partial \theta} \left( -\left(1 + \frac{1}{\rho^2}\right)\sin(\theta) - \frac{1}{\rho} \right) + 0 \right\}= \frac{1}{\rho} \left\{ (1 + \frac{1}{\rho^2})\cos(\theta) -\left(1 + \frac{1}{\rho^2}\right)\cos(\theta) \right\}= 0 [/math]

Como se puede observar también se demuestra que la divergencia es nula, dado que [math] \boxed {\nabla \cdot \vec u = 0} [/math] , lo que indica que el volumen del fluido es constante, es decir el fluido no se expande ni se contrae.

6 .-Líneas de corriente del Campo de Velocidades

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

[math]\vec{v}[/math] = [math]\vec{k}[/math] x [math]\vec{u}[/math] = [math]\begin{vmatrix} e_\rho & e_\theta & e_z \\ 0 & 0 & 1\\ cos(\theta)\cdot (1-\frac{1}{\rho ^2}) & sin(\theta) \cdot (1-\frac{1}{\rho^2}) -\frac{1}{\rho} & 0 \end{vmatrix}=[/math] [math] [(1+\frac{1}{\rho ^2}) \cdot sin(\theta) + \frac{1}{\rho}]e_\rho + [(1-\frac{1}{\rho ^2}) \cdot cos(\theta)]e_\theta[/math].


Como se ha comprobado en el apartado 5, el rotacional de [math]\vec{u}[/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 [/math]

[math]\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]\psi[/math] = [math] (\rho - \frac{1}{\rho}) \cdot sin(\theta) + Ln(\rho) [/math].
Una vez se obtiene tanto [math]\psi[/math] como [math]\vec{u}[/math], se introducen en MATLAB, así comprobando que efectivamente los vectores que componen [math]\vec{v}[/math] son ortogonales a las línes de corriente de [math]\psi[/math]. Para ello se ha empleado el siguiente código:

Lineas 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)((sin(theta).*(rho-1./rho))+log(rho));
Z=C(rho,theta);
contour(X,Y,Z,30);
hold on
%Grandiente de PSI
DX=((1+(1./(rho.^2))).*cos(theta).*sin(theta))+(cos(theta)./rho)+((-1+(1./(rho.^2))).*cos(theta).*sin(theta));
DY=((sin(theta).^2).*(1+(1./(rho.^2))))+(sin(theta)./rho)+((cos(theta).^2).*(1+(1./(rho.^2))));
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:

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);
C=@(x,y)((sin(theta).*(rho-1./rho))+log(rho));
Z=C(I,J);
contour(X,Y,Z,30);
hold on
X=rho.*cos(theta);      
Y=rho.*sin(theta);
%Grandiente de PSI
DX=((1+(1./(rho.^2))).*cos(theta).*sin(theta))+(cos(theta)./rho)+((-1+(1./(rho.^2))).*cos(theta).*sin(theta));
DY=((sin(theta).^2).*(1+(1./(rho.^2))))+(sin(theta)./rho)+((cos(theta).^2).*(1+(1./(rho.^2))));
quiver(X,Y,DX,DY);
%Gradiente de la función potencial
DXX=((cos(theta).^2).*(1-1./(rho.^2)))+((sin(theta).^2).*(1+ 1./(rho.^2)))+(sin(theta)./rho);
DYY=((sin(theta).*cos(theta)).*(1-1./(rho.^2)))+((sin(theta).*cos(theta)).*(-1-1./(rho.^2)))-(cos(theta)./rho);
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

Para una mayor apreciación, de las tangencias que forma [math] \vec u [/math] a las líneas de corriente y de la ortogonalidad formada por [math] \vec v [/math] también con las las líneas de corriente, se ha ampliado la fotografía de la gráfica anterior en un punto cualquiera, dado que se cumple a lo largo de todo el campo. Como se comprueba en la siguiente fotografía:

centro

7 .-Puntos de Remanso

En el borde del obstáculo [math]\rho = 1[/math] lo que quiere decir que [math]\vec{u}(1,\theta) = (-2 \cdot sin(\theta) - 1)[/math].
El módulo de [math]\vec{u}(1,\theta)[/math] es |[math]\vec{u}[/math]| = [math]\sqrt(4 \cdot sin(\theta)^2 + 4 \cdot sin(\theta) + 1)[/math]
Para que |[math]\vec{u}[/math]| alcance su valor máximo, [math]sin(\theta)[/math] tiene que ser igual a 1. Esto se alcanza en el ángulo [math]\frac{\pi}{2}[/math].
Para calcular el punto de remanso se ha de igualar |[math]\vec{u}[/math]| a 0 y despejar el valor de [math]\theta[/math]:

[math]\sqrt(4 \cdot sin(\theta)^2 + 4 \cdot sin(\theta) + 1) = 0[/math]

Finalmente [math]\theta = \frac{11 \pi}{6}[/math].

Para observar estos cálculos visualmente, una vez más se a codificado todo en MATLAB para así observar como en el punto de remanso efectivamente no se aprecia ningún vector:

Punto de Remanso

clear;clc;
I=linspace(1,5,50);   
J=linspace(0,2*pi,50);
[rho,theta]=meshgrid(I,J);  
X=rho.*cos(theta);      
Y=rho.*sin(theta);
f=@(x,y)(x+(1./x)).*cos(y)-y;
contour(X,Y,f(rho,theta),50);
hold on
DX=((cos(theta).^2).*(1-1./(rho.^2)))+((sin(theta).^2).*(1+ 1./(rho.^2)))+(sin(theta)./rho);
DY=((sin(theta).*cos(theta)).*(1-1./(rho.^2)))+((sin(theta).*cos(theta)).*(-1-1./(rho.^2)))-(cos(theta)./rho);
quiver(X,Y,DX,DY);
plot(1*cos(J),1*sin(J),'k','lineWidth',1);
%Punto de remanso previamente calculado
plot(0.866025,-0.5,'r*','LineWidth',3);
axis([-5,5,-5,5]);
colorbar;
title('Punto de Remanso');
xlabel('X');
ylabel('Y');
axis equal
view(2);
hold off

Para una mayor claridad se ha ampliado la gráfica anterior, para así conseguir una mayor nitidez del punto de remanso:

centro