Fluido Alrededor De Un Obstáculo Circular (Grupo 7)

Trabajo realizado por estudiantes
Título	Fluido alrededor de un obstáculo circular. Grupo 7
Asignatura	Teoría de Campos
Curso	2023-24
Autores	Sara Jorge García Álvaro Capilla Sanz Victor Del Águila Alonso Paula Gil Rodríguez Javier Moya López
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

En el presente artículo se analizará el comportamiento del flujo de un fluido incompresible alrededor de un obstáculo de forma circular, visualizando sus campos vectoriales y escalares. En el supuesto dado, el obstáculo se tratará de un anillo de radio 1, cuyo centro se encontrará en el origen de coordenadas y cuyo flujo estará definido por la siguiente función potencial:

[math]\varphi=(\rho+\frac{1}{\rho})cos(\theta)+\sqrt{2}\ln{\rho}[/math].

El fluido que rodea la superficie definida estará delimitado por: [math](\rho,\theta) \in [1,5]*[0,2\pi][/math].

2 Mallado

Para representar los puntos de la región ocupada por el fluido, se comenzará generando un mallado que describa el anillo situado entre los radios 1 y 5, con el origen de coordenadas como centro. Este proceso implica la creación de un conjunto de puntos distribuidos en coordenadas polares, el cual permitirá la visualización del mallado buscado. Se dibujarán los ejes en el intervalo [math](x,y) ∈ [-5,5]\times[-5,5] [/math].

Implementando el siguiente código en Matlab, se obtiene la superficie de trabajo.

rho=linspace(1,5,40);                        %definición de parámetros 
th=linspace(0,2*pi,40);

[U,V]=meshgrid(rho,th);                      %creación del mallado

hold on
X=U.*cos(V);                                 %parametrización de la superficie
Y=U.*sin(V);
Z=0.*U;
mesh(X,Y,0*Z);                               %representación de la placa 
plot(1*cos(th),1*sin(th),'k','lineWidth',2); %dibujo del obstáculo
axis([-5,5,-5,5]);                           %delimitación de ejes 
view(2);                                     %vista en proyección bidimensional
title ('PLACA');
xlabel 'ABSCISAS'
ylabel 'ORDENADAS'
axis equal
hold off

Representación del Mallado

3 Función Potencial y Campo de Velocidades

Tras definir la región de estudio, se analizará tanto la función potencial como el campo de velocidades que definen el comportamiento del fluido dado:

3.1 Función Potencial

Se procede a representar visualmente la función potencial del fluido, siendo esta:

[math]\varphi=(\rho+\frac{1}{\rho})cos(\theta)+\sqrt{2}\ln{\rho}[/math]

u=linspace(1,5,50);
v=linspace(0,2*pi,50);
[rho,th]=meshgrid(u,v);
Mx=rho.*cos(th);
My=rho.*sin(th);
Mz=zeros(size(Mx));
f=(rho+1./rho).*cos(th)+sqrt(2).*log(rho);
contour(Mx,My,f,25); colorbar
axis([-5,5,-5,5]); axis equal

3.2 Campo de Velocidades

La velocidad de un campo escalar se describe mediante el gradiente de dicha función. Así, el campo de velocidades se definirá como [math]\vec u=\nabla\varphi[/math], definiéndose el gradiente como:

[math]\vec u=\nabla\varphi=\frac{∂φ}{∂x}\vec i + \frac{∂φ}{∂j}\vec j [/math]

Para el cálculo del gradiente en coordenadas cartesianas, se tendrán en cuenta las siguientes igualdades matemáticas:

[math]\rho=\sqrt{x^2+y^2}\qquad\theta=arctg(\frac{y}{x})[/math]

Teniendo en cuenta todo lo descrito anteriormente, se tiene que el campo de velocidades de la función potencial dada será:

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

Habiendo calculado ya el campo de velocidades, se procede a hacer una representación gráfica de este mediante su visualización en el programa Matlab:

u=linspace(1,5,50);
v=linspace(0,2*pi,50);
[rho,th]=meshgrid(u,v);
Mx=rho.*cos(th);
My=rho.*sin(th);
Mz=zeros(size(Mx));
fx=@(x,y) -(1+((y.^2-x.^2)./((x.^2+y.^2).^2))+((sqrt(2).*x)./(x.^2+y.^2)));
fy=@(x,y) -(((-2)*(x.*y)./(x.^2+y.^2).^2)+((sqrt(2).*y)./(x.^2+y.^2)));   
colorbar;
quiver(Mx,My,fx(Mx,My),fy(Mx,My))
axis([-5,5,-5,5]); axis equal

4 Ortogonalidad entre el Gradiente ([math]\vec u[/math]) y un Vector Normal Asociado ([math]\vec n[/math])

En este caso [math]\vec{e_z}[/math] es un posible vector normal al campo de velocidades, ya que la función potencial dada existe en el plano [math]z=0[/math], por lo que el producto escalar será nulo.

En la frontera del obstáculo, [math]\rho=1 [/math], por lo que el producto escalar [math]\vec u\cdot\vec n=0[/math]

[math](\sqrt{2}\vec{e_\rho}-2sin(\theta)\vec{e_\theta})\cdot(\vec{e_z})=0[/math], siendo este primer vector el campo de velocidades del fluido particularizado en [math]\rho=1 [/math]

En la frontera del obstáculo, la velocidad es nula porque el vector velocidad será perpendicular al vector normal de la superficie. Por ello, las partículas de fluido que chocan con el obstáculo se detendrán de forma que rodearán el obstáculo para seguir avanzando.

5 Valor del Campo de Velocidades en Punto Lejano al Obstáculo

Sea [math]\rho[/math] un valor considerablemente elevado, se puede despreciar el término [math]\frac{1}{\rho}[/math], por lo que la expresión del campo de velocidades será la siguiente:

[math]\vec{u}=(cos(\theta)+\frac{\sqrt{2}}{\rho})\vec{e_\rho}-sin(\theta)\vec{e_\theta}[/math]

Acto seguido, se procede al cambio de coordenadas cilíndricas a cartesianas haciendo uso de la matriz cambio de base, como se muestra a continuación:

[math]\begin{pmatrix} u_1\\ u_2\\ u_3 \end{pmatrix}=\begin{pmatrix} \cos \theta & -sen \theta & 0\\ sen \theta & \cos \theta & 0\\ 0 & 0 & 1 \end{pmatrix}·\begin{pmatrix} \cos (\theta)+frac{\sqrt{2}}{\rho}\\ -sen \theta \\ 0 \end{pmatrix}=\begin{pmatrix} (cos(\theta)+\frac{\sqrt{2}}{\rho})cos(\theta)+sen^2(\theta)\\ (cos(\theta)+\frac{\sqrt{2}}{\rho})sin(\theta)-sen(\theta)cos(\theta)\\ 0 \end{pmatrix}=\begin{pmatrix} \frac{\sqrt{2}}{x^2+y^2}x\\ \frac{\sqrt{2}}{x^2+y^2}y\\ 0 \end{pmatrix}[/math]

6 Incompresibilidad

En el presente apartado se procederá a comprobar que tanto el rotacional como la divergencia del campo de velocidades son nulos, es decir, comprobar que se trata de un fluido incompresible.

6.1 Divergencia

La divergencia de un campo vectorial responderá a la siguiente fórmula:

Sea [math]\vec(F)=F_\rho\vec(e_\rho)+F_\theta\vec(e_\theta)+F_z\vec(e_z)[/math]

[math]\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}(\rho)(F_\rho)+\frac{\partial}{\partial{\theta}}(F_\theta)+\frac{\partial}{\partial{z}}(\rho)(F_z)][/math]

[math]\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}[(\rho-\frac{1}{\rho})cos(\theta)+\sqrt{2}]+\frac{\partial}{\partial{\theta}}[-(1+\frac{1}{\rho^2})sin(\theta)]]=\frac{1}{\rho}[cos(\theta)(1+\frac{1}{\rho^2})]+[-cos(\theta)(1+\frac{1}{\rho^2})]=\vec {0}[/math]

Que la divergencia resulte nula significa que el volumen del fluido se mantendrá constante (no se expandirá ni se contraerá)

6.2 Rotacional

La divergencia de un campo vectorial responderá a la siguiente fórmula:

Sea [math]\vec(F)=F_\rho\vec(e_\rho)+F_\theta\vec(e_\theta)+F_z\vec(e_z)[/math]

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

Que el rotacional resulte nulo significa que las partículas de fluido no giran

7 Líneas de Corriente

Se comenzará calculando el campo [math]\vec v[/math], que en cada punto es ortogonal a [math]\vec u[/math] ([math]\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)+\frac{\sqrt{2}}{\rho} & (1+\frac{1}{\rho^2})sin(\theta) & {0} \end{vmatrix}= (1+\frac{1}{\rho^2})sin(\theta)\vec {e}_{\rho} + [(1-\frac{1}{\rho^2})cos(\theta)+\frac{\sqrt{2}}{\rho}]\vec {e}_{\rho} =\vec v[/math]

A continuación, se comprobará que [math]\vec v[/math] es irrotacional. Así:

[math]\vec v\times\nabla= \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 & {0} \end{vmatrix}=\frac{1}{\rho}[[(1+\frac{1}{\rho^2})cos(\theta)]\vec {e}_{z}-[(1+\frac{1}{\rho^2})cos(\theta)]\vec {e}_{z}]=\vec {0}[/math]

Se procede a calcular [math]\psi[/math] resolviendo el sistema de ecuaciones [math]\nabla\cdot\psi=\vec v[/math]

[math]\frac{\partial\psi}{\partial\rho}=v_\rho=\int (1+\frac{1}{\rho^2})sen(\theta)\,d\rho=sen(\theta) (\rho-\frac{1}{\rho})+f(\theta)[/math]

[math]\frac{\partial\psi}{\partial\theta}=\rho v_\theta=\int (\rho-\frac{1}{\rho})cos(\theta)+sqrt{2}\,d\theta=sen(\theta) (\rho-\frac{1}{\rho})+\sqrt{2}\theta[/math]

De ello se obtendrá la siguiente igualdad, que representa el potencial escalar:

[math]\psi=sen(\theta) (\rho-\frac{1}{\rho})+\sqrt{2} (\theta)[/math]

A continuación, se representa el campo de velocidades y el potencial escalar mediante el siguiente código de Matlab

u=linspace(1,5,250);
v=linspace(0,2*pi,250);
[rho,th]=meshgrid(u,v);
Mx=rho.*cos(th);
My=rho.*sin(th);
Mz=zeros(size(Mx));


%[math]\psi[/math]
f=@(x,y) ((y)-(y./(x.^2+y.^2)))+(atan(y./x).*((2)^(1/2)));
%esto es u
f1=@(x,y) -(1+((y.^2-x.^2)./((x.^2+y.^2).^2))+((sqrt(2).*x)./(x.^2+y.^2)));
f2=@(x,y) -(((-2)*(x.*y)./(x.^2+y.^2).^2)+((sqrt(2).*y)./(x.^2+y.^2)));
hold on
subplot(1,1,1);
contour(Mx,My,f(Mx,My),100);
quiver(Mx,My,f1(Mx,My),f2(Mx,My));
axis equal
hold off

8 Puntos de Remanso

Se define punto de remanso como aquel lugar donde el módulo de la velocidad del flujo en nulo. Se procede a determinar si existe o no algún punto de remanso en el flujo dado:

9 Presión del Fluido

[math][/math]