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

Representación de la Función Potencial

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

Representación del Campo de Velocidades

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]

Se verifica representando gráficamente que [math]\vec v[/math] y el campo de velocidades [math]\vec u[/math] son ortogonales, adjuntando también el código de Matlab correspondiente:

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


%esto es v
fx=@(x,y) (((2).*(x.*y))./((x.^2+y.^2).^(2)))-((sqrt(2).*y)./((x.^2)+(y.^2)));
fy=@(x,y) 1+((sqrt(2).*x)./(x.^2+y.^2))+((y.^2-x.^2)./(x.^2+y.^2).^(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);
quiver(Mx,My,f1(Mx,My),f2(Mx,My));
quiver(Mx,My,fx(Mx,My),fy(Mx,My));
axis equal
hold off

Representación de las Líneas de Corriente y de las Curvas Equipotenciales

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

Representación de las Líneas de Corriente y de las Curvas Equipotenciales

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:

[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]

Particularizando en [math]\rho=1[/math], se tiene que el campo de velocidades valdrá:

[math]\vec u=\sqrt{2}\vec{e}_\rho-2sen(\theta)\vec{e}_\theta [/math]

Hallando ahora su módulo, se tiene que:

[math]|\vec u |=\sqrt{2+4sen^2(\theta)}[/math]

Se alcanzará el máximo cuando [math]sen(\theta)=1[/math], es decir, cuando [math]\theta=\frac{\pi}{2}[/math] ó [math]\theta=\frac{3\pi}{2} [/math] y el mínimo cuando [math]sen(\theta)=0[/math], es decir, cuando [math]\theta=0[/math] ó [math]\theta=\pi [/math]

Los puntos de remanso se hallan calculando:

[math]|\vec u|=0[/math]

Es decir:

[math]0=2+4sen^2(\theta)[/math]

Sin embargo, la igualdad no se cumple en ningún caso. Por lo tanto, se determina que no existen puntos de remanso.

9 Presión del Fluido

El principio de Bernoulli, en dinámica de fluidos, describe le comportamiento de una fluido moviéndose a lo largo de una línea de corriente. En la ecuación, se relacionan las siguientes variables:

[math]p :[/math]presión estática a la que está sometido el fluido debida a las moléculas que lo rodean.
[math]d :[/math]densidad del fluido.
[math]\vec {u} :[/math]velocidad de flujo del fluido.

El teorema de Bernoulli es el siguiente:

[math]\frac { 1 }{ 2 } d { \left| \vec { u } \right| }^{ 2 }+ { p } =cte[/math]

Para que esta expresión pueda aplicarse, el fluido debe ser incompresible, no viscoso y fluir en régimen estacionario (la velocidad en un punto determinado no varía con el tiempo). Se supondrá que el fluido cumple estas condiciones.
Suponiendo que el fluido efectivamente cumple la ecuación de Bernouilli, que posee una densidad constante de d=2 y tomando como cte=10, se calculará la presión del fluido.

Por lo tanto:

[math]p=10-{ \left| \vec { u } \right| }^{ 2 }[/math]

Calculamos el módulo al cuadrado de nuestro campo de velocidades:

[math]{ \left| \vec { u } \right| }=cos^2(\theta)(1+\frac{1}{\rho^4}-\frac{2}{\rho^2})+\frac{2}{\rho^2}+\frac{2\sqrt2}{\rho}cos(\theta)(1+\frac{1}{\rho^2})+sen^2(\theta)(1+\frac{1}{\rho^4}+\frac{2}{\rho^2})[/math]

Por lo tanto, la presión que define la presión del fluido es:

[math]p=10-[cos^2(\theta)(1+\frac{1}{\rho^4}-\frac{2}{\rho^2})+\frac{2}{\rho^2}+\frac{2\sqrt2}{\rho}cos(\theta)(1+\frac{1}{\rho^2})+sen^2(\theta)(1+\frac{1}{\rho^4}+\frac{2}{\rho^2})][/math]

[math][/math]

11 Trayectoria de la línea de corriente

Si se examina una partícula del fluido tratado siguiendo la trayectoria de una línea de corriente, se podría observar que el fluido comienza su trayectoria original y en cuanto se acerca al obstáculo, comienza a rodear su contorno. Una vez pasado el obstáculo, el fluido continuará con su trayectoria original. Si se relaciona la variación de velocidad y presión se comprueba, observando las gráficas, que cuando las presiones son máximas, las velocidades son mínimas. Se confirma que en los puntos donde la velocidad es 0 (calculado anteriormente), [math]\theta=0[/math] [math]\theta=\pi[/math] , la presión es máxima. Es decir, la presión y la velocidad de las partículas son inversamente proporcionales.