Flujo alrededor de un obstáculo circular (Grupo 37)

Trabajo realizado por estudiantes
Título	Flujo alrededor de un obstáculo circular. Grupo 37
Asignatura	Teoría de Campos
Curso	2025-26
Autores	Paula Gutiérrez Pascual Rafael Martín Candilejo Jaime Mateos Bermejo Hugo Zamora Ramos
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

El flujo es la cantidad de fluido que atraviesa una superficie por unidad de tiempo. Es decir, el flujo de un fluido nos marca el movimiento de este desde un lugar a otro.

Dadas las descripciones anteriores, no es difícil caer en la cuenta de que el flujo de un fluido será capaz de describirnos como este se desplaza a través de una sección de interés, siendo capaces de analizar la velocidad y dirección de su movimiento del fluido en cada punto; es decir, el campo de velocidades del mismo.

Si el interés se dirige a la mecánica de fluidos, podremos sacar jugosa información sobre los efectos internos al fluido estudiado, como la divergencia y rotacional. La divergencia nos marca el cambio del volumen del fluido al enfrentarse al movimiento, marcado por el campo vectorial de velocidades, mientras que el rotacional muestra la tendencia del fluido, de su campo vectorial, a rotar alrededor de un punto.

Aplicado al problema planteado, al tratarse de un fluido incompresible, el volumen siempre se conserva y, por consiguiente, la divergencia del mismo siempre será cero.

Otra información destacable que podemos sacar del campo de velocidades del fluido será la interacción con las paredes u obstáculos de la sección de interés. Esto tendrá una cierta trascendencia para el desarrollo de nuestro trabajo, dado que se plantea una situación donde el fluido estudiado interacciona con un obstáculo de forma circular.

Frente a esta situación, el desarrollo del trabajo se hará respecto a coordenadas cilíndricas.

1 Mallado

Se considera que el obstáculo mencionado coincide con el círculo unidad, con centro en el origen de coordenadas; luego, la región ocupada por el fluido será el exterior del círculo considerado.

El primer objetivo del grupo es la representación de los puntos interiores de la región ocupada por el fluido, para ello, con la ayuda del programa MATLAB, se dibujará un mallado que cumpla la representación.

Este mallado se ajusta al anillo de radio interior 1 y radio exterior 5 y, al igual que el obstáculo, centro en el origen de coordenadas.

Como ultima observación, para ilustrar que el fluido ocupa el exterior del obstáculo, los ejes se mostrarán en el intervalo [−4,4]×[−4,4].

Mallado

r=linspace(1,5,50);
a=linspace(0,2*pi,50);

[R,A]=meshgrid(r,a);

hold on
X=R.*cos(A);
Y=R.*sin(A);
Z=zeros(size(A));

mesh(X,Y,Z);

plot3(cos(A), sin(A), zeros(size(A)), 'k', 'LineWidth', 1);
hold off

axis equal;
axis([-4,4,-4,4]);

xlabel ('Eje X');
ylabel ('Eje Y');
title ('Mallado de la región del fluido');

La representación del mallado facilita el análisis e intuición del comportamiento del fluido. Esto se debe gracias a la división del espacio que ocupa en pequeñas celdas que se comportan como unidades de cálculo más manejables, lo que permite evaluar de manera precisa velocidades, temperaturas y otros fenómenos.

2 Velocidad del fluido

Sea la función potencial

[math] \varphi (\rho ,\theta, z)=(\rho +\frac{1}{\rho})\cos (\theta ) [/math]

Hallamos su función gradiente tal que [math]\vec{u}[/math]=∇φ.

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

Función potencial

r=linspace(1,5,50);
a=linspace(0,2*pi,50);

[R,A]=meshgrid(r,a);

hold on
X=R.*cos(A);
Y=R.*sin(A);
Z=(R+1./R).*cos(A);

surf(X,Y,Z);

plot3(cos(A), sin(A), zeros(size(A)), 'k', 'LineWidth', 2);
hold off

axis equal;
axis([-4,4,-4,4]);

xlabel('Eje X');
ylabel('Eje Y');
title('Función potencial');
colorbar;

Campo de velocidades y función potencial

Detalle del campo de velocidades

r=linspace(1,5,40);
a=linspace(0,2*pi,40);

[R,A]=meshgrid(r,a);

hold on
X1=R.*cos(A);
Y1=R.*sin(A);
Z1=(R+1./R).*cos(A);

contour(X1,Y1,Z1,50);

x2=(cos(A))-(cos(A)./R.^2);
y2=-(sin(A))-(sin(A)./R.^2);

X2=cos(A).*x2-(sin(A)./R).*y2;
Y2=sin(A).*x2+(cos(A)./R).*y2;

quiver(X1,Y1,X2,Y2,'m');

plot3(cos(A), sin(A), zeros(size(A)), 'k', 'LineWidth', 1);
hold off

axis equal;
axis([-4,4,-4,4]);

xlabel('Eje X');
ylabel('Eje Y');
title('Campo de velocidades');
colorbar;

3 Divergencia y rotacional

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

3.1 Rotacional nulo

[math] \nabla\times\vec{u} =\frac{1}{\rho} \begin{vmatrix} \vec{e}_\rho & \rho\vec{e}_\theta & \vec{e}_z \\ \dfrac{\partial}{\partial\rho} & \dfrac{\partial}{\partial\theta} & \dfrac{\partial}{\partial z} \\ \left(1-\dfrac{1}{\rho^2}\right)\cos\theta & -\left(1+\dfrac{1}{\rho^2}\right)\sin\theta & 0 \end{vmatrix} =(-1 + 1/\rho^{2})\sin\theta\,\vec{e}_{z} \;-\; (-1 + 1/\rho^{2})\sin\theta\,\vec{e}_{z} = 0 [/math]

3.2 . Comprobación de la divergencia nula

Conociendo la fórmula de la divergencia calculamos:

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

[math] \nabla\cdot\vec{u} = \frac{1}{\rho} \left[ \frac{\partial}{\partial\rho} \Bigl(\rho\, \left(1-\frac{1}{\rho^{2}}\right)\cos\theta \;\vec{e}_{\rho} \Bigr) \;-\; \frac{\partial}{\partial\theta} \Bigl( \left(1+\frac{1}{\rho^{2}}\right)\sin\theta \; \vec{e}_{\theta} \Bigr) \right]=\frac{1}{\rho} \left[ \frac{\partial}{\partial\rho} \Bigl(\left(\rho-\frac{1}{\rho}\right)\cos\theta \;\vec{e}_{\rho} \Bigr) \;-\; \frac{\partial}{\partial\theta} \Bigl( \left(1+\frac{1}{\rho^{2}}\right)\sin\theta \; \vec{e}_{\theta} \Bigr) \right] [/math]

[math] \nabla\cdot\vec{u} = \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]

4 Líneas de corriente

Primero calcularemos el campo [math]\vec{v}[/math], que en cada punto es ortogonal a [math]\vec{u}[/math], ([math]\vec{v}[/math] = [math]\vec{k}\times\vec{u}[/math], donde [math]\vec{k}[/math]=[math]\vec {e}_{z}[/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+\frac{1}{\rho^2})sin(\theta)\vec {e}_{\rho} + [(1-\frac{1}{\rho^2})cos(\theta)]\vec {e}_{\theta} =\vec v[/math]

Comprobamos que [math]\vec{v}[/math] es irrotacional:

[math]\nabla\times\vec v= \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}} \\ v_\rho & \rho v_\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]

A continuación calculamos [math]\psi[/math], para ello resolveremos 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})sin(\theta)\,d\rho=sin(\theta) (\rho-\frac{1}{\rho})+f(\theta)[/math]

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

[math]\psi = \sin(\theta)\left(\rho - \frac{1}{\rho}\right)[/math]

Líneas de corriente del campo de velocidades

r=linspace(1,5,20);
a=linspace(0,2*pi,20);

[R,A]=meshgrid(r,a);

hold on
X3=R.*cos(A);
Y3=R.*sin(A);
Z3=sin(A).*(R-(1./R));

contour(X3,Y3,Z3,50);

x2=(cos(A))-(cos(A)./R.^2);
y2=-(sin(A))-(sin(A)./R.^2);

X2=cos(A).*x2-(sin(A)./R).*y2;
Y2=sin(A).*x2+(cos(A)./R).*y2;

quiver(X3,Y3,X2,Y2,'m');

plot3(cos(A), sin(A), zeros(size(A)), 'k', 'LineWidth', 1);
hold off

axis equal;
axis([-4,4,-4,4]);

xlabel('Eje X');
ylabel('Eje Y');
title('Líneas de corriente del campo de velocidades');
colorbar;

5 Velocidades en la frontera de S

Dada nuestra función de velocidades del fluido u ya calculada anteriormente, calcularemos la velocidad máxima, mínima y nula en la frontera del obstáculo circular S el cual viene descrito por la circunferencia unidad centrada en (0,0) Puesto que nuestra función esta en coordenadas cilíndricas nos será mas fácil analizar la frontera S puesto que solo tendremos que sustituir rho por 1

• De esta manera [math]u(\theta) = -2\sin\theta[/math].

Al estar analizando los puntos donde la velocidad es máxima, mínima y los puntos de remanso estudiaremos el modulo de

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

Como nuestra función es el modulo del seno es fácil analizar los valores

5.1 Puntos de Velocidad Máxima :

Se dan cuando [math]|\sin\theta| = 1[/math].

• [math]\theta = \pi/2 [/math]

• [math]\theta = 3\pi/2 [/math]

5.2 Puntos de Remanso (Velocidad mínima = 0):

Se dan cuando [math]\sin\theta = 0[/math]. Adeemas por como es la funcion modulo de u los puntos de remanso coincidiran con los minimos

• [math]\theta = 0 [/math]

• [math]\theta = \pi [/math]

En estos puntos, el fluido choca contra el obstáculo y se detiene momentáneamente.

6 Presión del fluido

Partimos de la ecuación de Bernoulli para un fluido con densidad [math]\rho = 2[/math]:

[math]\frac{1}{2}\rho|\vec{u}|^2 + p = \text{cte}[/math]

Y siendo u nuestra función del campo velocidades [math] \vec{u}= \left(1-\frac{1}{\rho^2}\right)\cos\theta\,\vec{e}_\rho - \left(1+\frac{1}{\rho^2}\right)\sin\theta\,\vec{e}_\theta [/math]

Entonces hallamos que [math]|\vec{u}|^2=1+\dfrac{1}{\rho^{4}}-\dfrac{2}{\rho^{2}}\cos2\theta[/math]

Sustituyendo la densidad por dos, utilizando una cte cualquira puesto que para el calculo de máximos y mínimos no nos es importante y despejando p llegamos a la ecuación

[math]p(\rho,\theta)=-\dfrac{1}{\rho^{4}}+\dfrac{2}{\rho^{2}}\cos2\theta[/math]

Para hallar los máximos y los mínimos de p derivaremos respecto a sus variables e igularemos a 0 para resolver el sistema, por último utilizaremos el teorema de weistrass para analizar los máximos y los mínimosen la frontera de S

6.1 Máximos y mínimos absolutos

[math] \begin{cases} \frac{\partial p}{\partial \rho} = 0 \\ \frac{\partial p}{\partial \theta} = 0 \end{cases} [/math]

[math] \begin{cases} \frac{4}{\rho^5} - \frac{4 \cos(2\theta)}{\rho^3}= 0 \\ \frac{4 \sin(2\theta)}{\rho^2}=0 \end{cases} [/math]

Resolviendo esto nos quedan cuatro soluciones para [math](\rho,\theta)[/math]

• [math](\rho,\theta)=(1,0)[/math]
• [math](\rho,\theta)=(1,\pi)[/math]

o Ahora para identificar si son máximos, mínimos o puntos de silla, utilizaremos la matriz hessiana

[math] H(f) = \begin{bmatrix} \frac{\partial^2 f}{\partial x_1^2} & \frac{\partial^2 f}{\partial x_1 \partial x_2} & \cdots & \frac{\partial^2 f}{\partial x_1 \partial x_n} \\ \frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial x_2^2} & \cdots & \frac{\partial^2 f}{\partial x_2 \partial x_n} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^2 f}{\partial x_n \partial x_1} & \frac{\partial^2 f}{\partial x_n \partial x_2} & \cdots & \frac{\partial^2 f}{\partial x_n^2} \end{bmatrix}[/math]

Haciendo las derivadas y sustituyendo quedaria

[math] H(p)(\rho,\theta)= \begin{pmatrix} \dfrac{12 \cos(2\theta)}{\rho^4} - \dfrac{20}{\rho^6} & \dfrac{8 \sin(2\theta)}{\rho^3} \\ \dfrac{8 \sin(2\theta)}{\rho^3} & -\dfrac{8 \cos(2\theta)}{\rho^2} \end{pmatrix}[/math] Ahora debemos analizar su determinante el cual es y sustituir por los valores encontrado para [math](\rho,\theta)[/math] [math]\det(H(p)) = \frac{-32 \rho^2 \cos^2(2\theta) - 64 \rho^2 + 160 \cos(2\theta)}{\rho^8}[/math]

Puesto que para [math](\rho,\theta)=(1,0) ,(1,\pi)[/math] el determinante queda positivo y la la primera componente de la primera columna y primera fila es negativo podemos concluir que [math](\rho,\theta)=(1,0),(1,\pi)[/math] son máximos absolutos

6.2 Máximos y mínimos en la frontera S

Dado que estamos trabajando en coordenadas cilíndricas y que la frontera S es realmente fácil de describir en estas coordenadas puesto que el radio no varia y solo varia el ángulo utilizaremos el teorema de Weistrass el cual establece que una función continua en un intervalo cerrado y acotado [a, b] siempre alcanza al menos un máximo absoluto y un mínimo absoluto dentro de ese intervalo. Esto significa que la función tiene un valor más alto y uno más bajo dentro de ese intervalo específico, garantizado por la continuidad en el conjunto cerrado y acotado, siendo nuestro intervalo [math]\rho \in [0, \pi], \quad \theta \in [1, \infty)[/math]

Por tanto sustituiremos por 1 el radio en la funcion de la presión [math]p(\rho,\theta)=-\dfrac{1}{\rho^{4}}+\dfrac{2}{\rho^{2}}\cos2\theta[/math] y nos quedara [math]p(\theta)=-1+2cos2\theta[/math]

Por tanto veremos que alcanza los máximos cuando [math]\theta=0,\pi[/math] y alcanza los máximos cuando [math]\theta=2\pi,2\pi/3[/math]

Máximos

• [math](\rho,\theta)=(1,0)[/math]
• [math](\rho,\theta)=(1,\pi)[/math]

Mínimos

• [math](\rho,\theta)=(1,\pi/2)[/math]
• [math](\rho,\theta)=(1,3\pi/2)[/math]

Observamos que los máximos en la frontera de S coindicen con los máximos absolutos calculados en el apartado anterior

6.3 Comparación con velocidad

Si nos remontamos al apartado donde calculamos las velocidades máximas y mínimas podemos observar que cuando la velocidad es máxima la presión es mínimas y a la inversa, esto se debe a la ecuación de Bernuilli que dicta que presión velocidad y altura guardan una relación llamada trinomio de Bernuilli el cual estudiaremos en profundida en el siguiente apartado

9 Apartado nueve

Sea la función potencial

[math] \varphi (\rho ,\theta, z)=(\rho +\frac{1}{\rho})\cos (\theta) +\frac{\theta}{4 \pi} [/math]

Función potencial

r=linspace(1,5,50);
a=linspace(0,2*pi,50);

[R,A]=meshgrid(r,a);

hold on
X=R.*cos(A);
Y=R.*sin(A);
Z=((R+1./R).*cos(A))+(A./4.*pi);

surf(X,Y,Z);

plot3(cos(A), sin(A), zeros(size(A)), 'k', 'LineWidth', 2);
hold off

axis equal;
axis([-4,4,-4,4]);

xlabel('Eje X');
ylabel('Eje Y');
title('Función potencial');
colorbar;

Hallamos su función gradiente tal que [math]\vec{u}[/math]=∇φ.

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

Campo de velocidades y función potencial

Detalle del campo de velocidades

r=linspace(1,5,40);
a=linspace(0,2*pi,40);

[R,A]=meshgrid(r,a);

hold on
X1=R.*cos(A);
Y1=R.*sin(A);
Z1=((R+1./R).*cos(A))+(A./4.*pi);

contour(X1,Y1,Z1,50);

x2=(cos(A))-(cos(A)./R.^2);
y2=-(sin(A))-(sin(A)./R.^2)+(1./4.*pi);

X2=cos(A).*x2-(sin(A)./R).*y2;
Y2=sin(A).*x2+(cos(A)./R).*y2;

quiver(X1,Y1,X2,Y2,'m');

plot3(cos(A), sin(A), zeros(size(A)), 'k', 'LineWidth', 1);
hold off

axis equal;
axis([-4,4,-4,4]);

xlabel('Eje X');
ylabel('Eje Y');
title('Campo de velocidades');
colorbar;

Con la fórmula del campo de velocidades,

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

calcularemos el rotacional y la divergencia.

9.1 Rotacional

[math] \nabla\times\vec{u} =\frac{1}{\rho} \begin{vmatrix} \vec{e}_\rho & \rho\vec{e}_\theta & \vec{e}_z \\ \dfrac{\partial}{\partial\rho} & \dfrac{\partial}{\partial\theta} & \dfrac{\partial}{\partial z} \\ \left(1-\dfrac{1}{\rho^2}\right)\cos\theta & -\left(1+\dfrac{1}{\rho^2}\right)\sin\theta & 0 \end{vmatrix} =(-1 + 1/\rho^{2})\sin\theta\,\vec{e}_{z} \;-\; (-1 + 1/\rho^{2})\sin\theta\,\vec{e}_{z} = 0 [/math]

9.2 Divergencia

Conociendo la fórmula de la divergencia calculamos:

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

[math] \nabla\cdot\vec{u} = \frac{1}{\rho} \left[ \frac{\partial}{\partial\rho} \Bigl(\rho\, \left(1-\frac{1}{\rho^{2}}\right)\cos\theta \;\vec{e}_{\rho} \Bigr) \;-\; \frac{\partial}{\partial\theta} \Bigl( \left(1+\frac{1}{\rho^{2}}\right)\sin\theta \; \vec{e}_{\theta} \Bigr) \right]=\frac{1}{\rho} \left[ \frac{\partial}{\partial\rho} \Bigl(\left(\rho-\frac{1}{\rho}\right)\cos\theta \;\vec{e}_{\rho} \Bigr) \;-\; \frac{\partial}{\partial\theta} \Bigl( \left(1+\frac{1}{\rho^{2}}\right)\sin\theta \; \vec{e}_{\theta} \Bigr) \right] [/math]

[math] \nabla\cdot\vec{u} = \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]

NO BORRAR