Fluido alrededor de un obstáculo circular (15-B)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Fluido alrededor de un obstáculo circular. Grupo 15-B |
| Asignatura | Teoría de Campos |
| Curso | 2022-23 |
| Autores | Julián Lavandero Pacheco Javier Sesmero Zamarrón Ana Sarró Redel Sergio Navarro Czyz Kevin Rosales Zambrana |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
1 Introducción
Un fluido es todo cuerpo que tiene la propiedad de fluir, carece de rigidez y elasticidad, y en consecuencia cede inmediatamente a cualquier fuerza tendente a alterar su forma, siendo capaz de adoptar la forma del recipiente que lo contiene. Los fluidos pueden ser líquidos o gases según la diferente intensidad de las fuerzas de cohesión existentes entre sus moléculas. En este artículo, se estudiarán los fluidos incompresibles y para ello, debemos conocer qué es un fluido incompresible. Estos, son aquellos en los que su volumen no disminuye al ejercerle fuerzas muy grandes, son, por tanto, los líquidos. Iniciaremos el estudio realizando un mallado.
2 Mallado
Vamos a comenzar dibujando un mallado que represente los puntos interiores de la región ocupada por un fluido. Este será el exterior del círculo unidad. Para la realización de este tomaremos un mallado del anillo comprendido entre los radios 1 y 5 y, centro, el origen. Para ilustrar que el fluido ocupa la parte exterior de un círculo, dibujaremos los ejes en el intervalo [math](x,y) ∈ [-5,5]\times[-5,5] [/math]
rho=linspace(1,5,30); %definimos rho y theta
th=linspace(0,2*pi,30);
[U,V]=meshgrid(rho,th); %Creamos la malla
hold on
X=U.*cos(V); %Parametrizamos la superficie
Y=U.*sin(V);
Z=0.*U;
mesh(X,Y,0*Z); %Dibujamos la placa
plot(1*cos(th),1*sin(th),'k','lineWidth',1); %Representación del obstáculo
axis([-5,5,-5,5]); %Establecemos los ejes
view(2);
title ('Placa');
xlabel 'EJE X'
ylabel 'EJE Y'
axis equal
hold off
3 Función potencial
Una vez definida la región de estudio, vamos a estudiar cómo se mueve el fluido y para ello, analizaremos la velocidad de las partículas de una función potencial dada. La velocidad de un campo escalar es el gradiente de dicha función. Para mayor claridad, vamos a definir brevemente un campo escalar y un campo vectorial.
Un campo escalar, definido en un dominio determinado [math] D \subset \mathbb{R}^{3} [/math] , es una aplicación que asigna a cada punto de D un escalar. Por otra parte, un campo vectorial definido en un dominio [math] D \subset \mathbb{R}^{3} [/math] es una aplicación que asigna a cada punto de D un vector libre.
Primero, realizaremos una representación de la función potencial del fluido, siendo esta [math] \varphi (\rho ,\theta)=(\rho +\frac{1}{\rho})\cos (\theta )+\sqrt{2}\theta ) [/math] . Posteriormente, representaremos el campo de velocidades.
rho=linspace(1,5,30); %Definimos rho y theta
th=linspace(0,2*pi,30);
[U,V]=meshgrid(rho,th); %Creamos la malla
hold on
X=U.*cos(V); %Parametrización de la superficie
Y=U.*sin(V);
f=@(rho,th)(rho+(1./rho)).*cos(th) +sqrt(2).*th; %Definimos la función potencial
Z=f(U,V); %Aplicamos la función potencial
surf(X,Y,Z); %Dibujamos la función
plot(1*cos(th),1*sin(th),'k','lineWidth',3); %Representamos el obstáculo
view(2);
axis([-5,5,-5,5]);
colorbar;
title ('Función potencial');
xlabel ('EJE X');
ylabel ('EJE Y');
axis equal
hold off
4 Campo de velocidades
La velocidad de dicha función potencial viene dada por su gradiente. Este representa el valor junto con la dirección de máximo crecimiento. Para su cálculo emplearemos la siguiente fórmula: [math] \vec u=\nabla \varphi=\left ( 1-\frac{1}{\rho^2} \right )\cdot\cos(\theta)\vec e_\rho - \frac{1}\rho\left [ \sin(\theta)\cdot(\rho+\frac{1}{\rho})+\sqrt 2 \right ]\vec e_\theta [/math]
rho=linspace(1,5,30); %Definimos rho y theta
th=linspace(0,2*pi,30);
[U,V]=meshgrid(rho,th); %Creamos la malla
X=U.*cos(V); %Parametrizamos la superficie
Y=U.*sin(V);
f=@(rho,th)(rho+(1./rho)).*cos(th) + sqrt(2).*th; %Definimos la función potencial
Z=f(U,V); %Aplicamos la función
contour(X,Y,Z,30); %Dibujamos las curvas de nivel
hold on
%Definimos las componentes X e Y del gradiente
Cx=(1-(1./U.^2)).*cos(V).^2 + (1+(1./U.^2)).*sin(V).^2 - (sqrt(2)./U).*sin(V);
Cy=(1-(1./U.^2)).*sin(V).*cos(V) - (1+(1./U.^2)).*sin(V).*cos(V) + (sqrt(2)./U).*cos(V);
quiver(X,Y,Cx,Cy); %Dibujamos el campo de velocidades
plot(1*cos(th),1*sin(th),'k','lineWidth',1); %Representación del obstáculo
axis([-5,5,-5,5]);
colorbar; %Añadimos una barra de color
axis equal
view(2);
hold off
4.1 Interpretación
Es importante observar mediante la ayuda de la gráfica, que, el campo de velocidades es ortogonal a las curvas de nivel de [math] \varphi [/math] . Una breve comprobación matemática de esta afirmación, sería considerar un vector n cualquiera normal a los puntos del obstáculo, y confirmando que, efectivamente, [math]\vec{u}\cdot \vec{n}=0 [/math]. Esto implica que son ortogonales entre sí.
Un vector ortogonal a las curvas de nivel sería [math]-\vec e_\rho [/math].
Evaluando el gradiente en [math]\rho = 1[/math], y posteriormente multiplicándolo por dicho vector perpendicular, podemos comprobar la condición de ortogonalidad:
[math]-\vec e_\rho \cdot \left ( -2\sin (\theta) +\sqrt{2} \right )\vec e_\theta = 0 [/math]
Si evaluamos [math] \vec u [/math] en la frontera del obstáculo obtendríamos lo siguiente:
[math]\vec u(1,\theta,z)=\nabla \varphi(1,\theta,z)= - 2\cdot\sin(\theta)+\sqrt 2 \vec e_\theta [/math]
4.2 Campo de velocidades lejos del obstáculo
Vamos a realizar una breve suposición para un mejor estudio de la función. Si estamos lejos del obstáculo, [math]\rho [/math] es muy grande, y podemos suponer que [math] \frac{1}{\rho } [/math] es despreciable. Si quisiésemos conocer el valor de u en esos puntos, operaríamos de la siguiente manera:
[math] \varphi = \rho \cos (\theta) +\sqrt{2}\theta [/math] , habiendo despreciado [math] \frac{1}{\rho } [/math].
Calculamos el gradiente de la función potencial resultante de despreciar [math] \frac{1}{\rho} [/math]
[math] \vec{u}= \nabla\varphi = \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= [/math]
[math]=\cos (\theta) \vec e_\rho +\left [ \frac{1}{\rho }\left ( -\sin (\theta) \cdot \rho \right )+\sqrt{2} \right ]\vec e_\theta =[/math]
[math] =\mathbf{\cos \vec e_\rho +\left ( \frac{\sqrt{2}}{\rho }-\sin (\theta) \right )\vec e_\theta} [/math]
4.3 Rotacional y divergencia
Procedemos a estudiar el rotacional y la divergencia del campo para conocer la naturaleza física del fluido.
El rotacional es un operador que mide la rotación en el movimiento de un fluido descrito por un campo vectorial de tres dimensiones. Muestra por tanto la tendencia de un campo a inducir rotación alrededor de un punto.
[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}\\ \cos \left (\theta \right ) \left ( 1-\frac{1}{\rho ^2} \right )& \rho \left [ -\frac{1}{\rho }sin \left ( \theta \right )\cdot\left ( \rho +\frac{1}{\rho } \right ) +\frac{\sqrt{2}}{\rho }\right ] & 0 \end{vmatrix}=[/math]
[math]=\frac{1}{\rho }\left [ \frac{\partial }{\partial \rho }\left [ -\sin \left ( \theta \right )\left ( \rho +\frac{1}{\rho } \right ) +\sqrt{2}\right ]\vec e_z -\frac{\partial }{\partial \theta }\left [ \cos \left ( \theta \right )\left ( 1-\frac{1}{\rho ^2 } \right ) \right ]\vec e_z\right ]=[/math]
[math]=\frac{1}{\rho }\left [ -\sin \left ( \theta \right ) \left ( 1-\frac{1}{\rho ^2} \right )\vec e_z+\sin \left ( \theta \right )\left ( 1-\frac{1}{\rho ^2} \right )\vec e_z\right ]=\mathbf{0} [/math]
Por el contrario, la divergencia de un campo vectorial es una cantidad escalar que mide la diferencia entre el flujo entrante y el flujo saliente en una superficie.
[math] \nabla\cdot \vec{u} =\frac{1}{\rho} \left [ \frac {\partial}{\partial \rho } \left ( \rho\cdot cos (\theta) \cdot(1 -\frac {1}{\rho^2}) \right ) + \frac {\partial}{\partial \theta} \left ( -\frac {1}{\rho}\cdot sin(\theta) \cdot (\rho + \frac {1}{\rho})+ \frac {\sqrt 2}{\rho} \right ) \right ] = [/math]
[math] =\frac{1}{\rho} \left [ \frac {\partial}{\partial \rho } \left ( cos (\theta) \cdot(\rho -\frac {1}{\rho}) \right ) + \frac {\partial}{\partial \theta} \left ( - sin(\theta) \cdot (1 + \frac {1}{\rho^2})+\frac {\sqrt 2}{\rho} \right ) \right ] =[/math]
[math] =\frac{1}{\rho} \left [ cos (\theta) \cdot (1 +\frac {1}{\rho^2}) - (1 + \frac {1}{\rho^2})\cdot cos(\theta) \right ] = \frac{1}{\rho}\cdot 0 = \textbf{0} [/math]
Obtenemos que el rotacional es nulo. Esto quiere decir que las partículas del fluido no giran. Del mismo modo, la divergencia es nula. Esto implica que el fluido mantiene su volumen sin producirse un intercambio de flujo. Hemos podido comprobar que, en efecto, estamos ante un fluido incompresible.
5 Líneas de corriente
Vamos a proceder a calcular las líneas de corriente del campo [math]\vec{u}[/math], es decir las líneas que son tangentes a [math]\vec{u}[/math] en cada punto. Estas líneas determinan la trayectoria que siguen las partículas del fluido.
Para ello calculamos el campo que en cada punto es ortogonal a [math]\vec{u}[/math]. Este vector lo llamaremos [math]\vec{v}[/math] y lo calcularemos de la forma [math]\vec{v}=\vec{k}\times \vec{u}[/math].
Como hemos visto en el apartado anterior, el rotacional al ser nulo implica que [math]\vec{v}[/math] sea irrotacional y además tiene un potencial escalar [math]\psi [/math] cuyo gradiente es [math]\vec{v}[/math]. Este potencial se conoce como función de corriente de [math]\vec{u}[/math].
Para calcular [math]\psi [/math] resolveremos el sistema de ecuaciones [math]\nabla\psi =\vec{v}[/math], y dibujaremos las líneas [math]\psi =cte[/math].
Comenzamos calculando el vector [math]\vec{v}[/math]
[math]\overrightarrow{v}=\overrightarrow{k}\times \overrightarrow{u}=\begin{vmatrix} \overrightarrow{e_\rho }& \overrightarrow{e_\theta } & \overrightarrow{e_z}\\ 0& 0& 1\\ \cos \left ( \theta \right )\left [ 1-\frac{1}{\rho ^2} \right ] & \frac{1}{\rho }\left [ -\sin \left ( \theta \right ) \left ( \rho +\frac{1}{\rho } \right )+\sqrt{2}\right ] & 0 \end{vmatrix}=[/math]
[math]=\left ( \frac{1}{\rho } \sin \left ( \theta \right ) \left ( \rho +\frac{1}{\rho } \right )-\frac{\sqrt{2}}{\rho }\right )\overrightarrow{e_\rho }+\cos \left ( \theta \right )\left ( 1-\frac{1}{\rho ^2} \right )\overrightarrow{e_\theta }[/math]
Resolvemos por tanto las siguientes ecuaciones:
[math]\frac{\partial \psi }{\partial \rho }=V_\rho [/math]
[math]\frac{1}{\rho }\frac{\partial \psi }{\partial \theta }=V_\theta [/math]
[math]\frac{\partial \psi }{\partial z }=V_z[/math]
Finalmente, obtenemos la línea de corriente cuyo valor es:
[math] \psi= sin (\theta)\cdot(\rho-\frac{1}{\rho}) - \sqrt 2 \cdot \ln \rho [/math]
Posteriormente comprobaremos que, efectivamente son las líneas de corriente de [math]\vec{u}[/math] siendo estas tangentes a la velocidad del fluido. Del mismo modo comprobamos que las líneas de corriente son ortogonales a las curvas equipotenciales.
rho=linspace(1,5,30);
th=linspace(0,2*pi,30);
[U,V]=meshgrid(rho,th); %Creamos la malla
X=U.*cos(V); %Parametrizamos la superficie
Y=U.*sin(V);
%dibujamos las lienas de corriente
psi=@(rho,th)(sin(V).*(U-(1./U)))- (sqrt(2).*log(U)); %Definimos la función ψ (psi)
Z=psi(U,V); %Aplicamos la función
contour(X,Y,Z,30); %Dibujamos las líneas de corriente
hold on
%Dibujamos el campo de velocidades (gradiente de phi)
Cx=(1-(1./U.^2)).*cos(V).^2 + (1+(1./U.^2)).*sin(V).^2 - (sqrt(2)./U).*sin(V);
Cy=(1-(1./U.^2)).*sin(V).*cos(V) - (1+(1./U.^2)).*sin(V).*cos(V) + (sqrt(2)./U).*cos(V);
quiver(X,Y,Cx,Cy);
plot(1*cos(th),1*sin(th),'k','lineWidth',1); %Representación del obstáculo
axis([-5,5,-5,5]);
colorbar;
xlabel 'Eje X';
ylabel 'Eje Y';
axis equal
view(2);
hold off
En la imagen de superior se puede observar que las líneas de corriente de [math] \vec u [/math] son tangentes a la velocidad del fluido.
En las imágenes inferiores se pueden observar las líneas de corriente de [math] \vec u [/math], las flechas naranjas representan el gradiente de [math] \varphi [/math], que a su vez es tangente a las líneas de corriente, y las flechas azules representan el gradiente de [math] \psi [/math], que a su vez es tangente a las curvas equipotenciales de [math] \varphi [/math].
Zoom de la gráfica superior
rho=linspace(1,5,30);
th=linspace(0,2*pi,30);
[U,V]=meshgrid(rho,th); %Creamos la malla
X=U.*cos(V); %Parametrizamos la superficie
Y=U.*sin(V);
hold on
psi=@(rho,th)(sin(V).*(U-(1./U)))- (sqrt(2).*log(U)); %Definimos la función ψ (psi)
Z=psi(U,V); %Aplicamos la función
contour(X,Y,Z,15); %Dibujamos las líneas de corriente
%dibujamos el campo de velocidades (gradiente de psi)
CX=(1+(1./U^2)).*cos(V).*sin(V) - (sqrt(2)./U).*cos(V) - (1-(1./U.^2)).*sin(V).*cos(V);
CY=(1+(1./U^2)).*sin(V).^2 -(sqrt(2)./U).*sin(V) + (1-(1./U.^2)).*cos(V).^2;
quiver(X,Y,CX,CY);
%Dibujamos el campo de velocidades ( gradiente de phi)
Cx=(1-(1./U.^2)).*cos(V).^2 + (1+(1./U.^2)).*sin(V).^2 - (sqrt(2)./U).*sin(V);
Cy=(1-(1./U.^2)).*sin(V).*cos(V) - (1+(1./U.^2)).*sin(V).*cos(V) + (sqrt(2)./U).*cos(V);
quiver(X,Y,Cx,Cy);
plot(1*cos(th),1*sin(th),'k','lineWidth',1); %Representación del obstáculo
axis([-5,5,-5,5]);
colorbar;
xlabel 'Eje X';
ylabel 'Eje Y';
axis equal
view(2);
hold off
6 Velocidades máxima y mínimas en un fluido
Calcular los puntos de la frontera del obstáculo S donde el módulo de la velocidad es mayor y menor. Los puntos en los que la velocidad es nula se denominan puntos de remanso. Señalar los puntos de remanso en el borde del obstáculo en una gráfica.
Puesto que [math] \vec{u} =\nabla\varphi [/math] tenemos que hallar el módulo de la velocidad, [math] \left \| \vec u \right \|[/math] . Esto es un escalar y puesto que se quiere analizar en la región S =[math] \left \{(\rho,θ): \rho=1 \right \}[/math] simplificando el módulo en la siguiente ecuación: [math] \left \| \vec u \right \|= -2\cdot sin (\theta) + \sqrt 2 [/math]
Derivando la ecuación se halla los puntos de inflexión obtiendose los máximos y mínimos. [math] \frac {\partial}{\partial\theta}\left ( -2sin\theta + \sqrt 2 \right )=0 \Rightarrow cos(\theta)=0 \Rightarrow \theta=3π/2 [/math] y [math]\theta=π/2 [/math]
Por lo consiguiente, los puntos de inflexión son: 3π/2 y π/2. Puesto que el módulo de la velocidad es dado en valor absoluto, su valor mínimo posible es [math] \left \| \vec u \right \|=0 [/math]
Dichos puntos son en 3π/4 y π/4
En conclusión: el máximo se da en: 3π/2 y el mínimo se da en: 3π/4 y π/4.
Observación: para comprobar que son máximos y mínimos se dan valores.
x=linspace(0,2*pi,1000);
y=f(x);
a=1;
%Se cambian de signo los valores negativos
for i=y
if i<0
i=-i;
end
y(a)=i;
a=a+1;
end
plot(x,y)
xlabel('ángulo')
ylabel('módulo de la velocidad')
7 Ecuación de Bernoulli
En esta sección, se calculará los puntos de mayor y menor presión del fluido. Supondremos que la densidad del fluido es constante [math] d = 2 [/math], y que se verifica la ecuación de Bernoulli: [math] 1/2 d \left \| \vec u \right \|{^2}=cte [/math]. La presión del fluido es calculada dando el valor 10 a la constante anterior.
Puesto que al valor de la cte se le asigna el valor de 10, la ecuación queda como: [math] p=10-\left \| \vec u \right \|{^2} [/math]. Dado que el interés se concentra en la región del obstáculo S, podemos dejar [math] \left \| \vec u \right \|{^2} [/math] en función del parámetro [math] \theta [/math].
[math] p=10-\left ( -2senθ +\sqrt[]{2} \right ){^2} [/math] Derivando la función e igualando a 0 se obtiene la siguiente ecuación: [math] 2senθcosθ=\sqrt[]{2}cosθ [/math]. Hallamos los distintos máximos y mínimos. Puesto que cosθ puede ser =0, debemos acotar la función, Esto significa que se tiene que comprobar que tanto en π/2 y 3π/2 pueda haber un máximo o mínimo.
El máximo se da en π/4 y 3π/4 mientras que el mínimo se da en: 3π/2.
Para observar las diferencias, a continuación se expone un gráfico de las presiones en el obstáculo.
tt=linspace(0,2*pi,100)
p=10-(2.*sin(tt)+sqrt(2)).^2
plot(tt,p)
Los puntos del obstáculo donde se alcanzan máximos de presión son opuestos a los puntos donde el módulo de la velocidad es mínimo y viceversa. En el gráfico se ha rotado media vuelta los valores de la presión para mejor visualización.
f=@(x) -2*sin(x) + sqrt(2);
x=linspace(0,2*pi,1000);
y=f(x);
a=1;
%Se cambian de signo los valores negativos
for i=y
if i<0
i=-i;
end
y(a)=i;
a=a+1;
end
hold on
plot(x,y);
tt=linspace(0,2*pi,100);
p=10-(2.*sin(tt+pi)+sqrt(2)).^2;
plot(tt,p)
xlabel('ángulo')
ylabel('módulo de la velocidad y presión')
hold off
8 Velocidades y presiones
Observando tanto el campo de velocidades como el de presiones, se puede deducir que hay una relación entre presión y velocidad. Esto se deduce a partir de que cuanto menor es su presión (parte baja del gráfico), aumenta el módulo de la velocidad. Además, esto es lógico ya que cuando se habla de un fluido incompresible, su valor de Trinomio de Bernoulli ha de ser cte.
B=valor del trinomio de Bernoulli
p= presión
z= altura
v=velocidad que lleva el fluido en determinado punto
g=valor de la gravedad
γ=valor del peso específico
[math] B=z+p/γ+v{^2}/(2g) [/math]
Luego, para compensar, si la presión desciende, su velocidad ha de aumentar y viceversa.
9 Paradoja de D'Alembert
Puesto que p [math] \vec {n} [/math] es la fuerza que ejerce el fluido en cada punto de la frontera, esto significa que, al sumar la proyección de todas estas fuerzas sobre la dirección [math] \vec {i} [/math] la resultante es nula. En particular, el fluido no realiza ningún empuje sobre el obstáculo en la dirección horizontal.
[math] \int p·\vec {n}·\vec {i} = \int p cosθ ds =\int 10-{(-2\cdot sin (\theta) + \sqrt 2)^2}cosθ\left \|γ′\right \| dθ[/math]
Si parametrizamos a γ como [math] \left \{ (\rho,\theta)=(1,\theta) \right \} \theta\epsilon[0,2π) \Rightarrow \left | γ′ \right |=1[/math]
Por lo que:
[math] \int_{0}^{2π} (10cosθ-2cosθ-4sen{^2}θcosθ+4\sqrt{2}senθcosθ)dθ [/math]
Al ser una integral, es lineal por lo que se puede descomponer en suma de varias integrales.
[math] \int_{0}^{2π} 10cosθdθ - \int_{0}^{2π} 2cosθdθ - \int_{0}^{2π} 4sen{^2}θcosθdθ + \int_{0}^{2π} 4\sqrt{2}senθcosθdθ [/math]
Todas provienen de primitivas del seno y al estar acotadas de 0 a 2π,
[math] \int_{0}^{2π} 10cosθdθ - \int_{0}^{2π} 2cosθdθ - \int_{0}^{2π} 4sen{^2}θcosθdθ + \int_{0}^{2π} 4\sqrt{2}senθcosθdθ = 0 [/math]
Esto es lo que se conoce como la paradoja de D’Alembert.
10 Curvas de nivel
Las curvas de nivel muestran un mínimo absoluto en el punto en coordenadas cilindricas (1, 3π/2, 0), y tres máximos en los puntos (1, π/4, 0), (1, 3π/4, 0) y en un punto por encima del obstáculo. Basándonos únicamente en las lineas de nivel resulta muy difícil hallar cuales de ellos son máximos absolutos y cuales máximos locales. Con un análisis más general, vemos que las presiones en la zona superior son mayores a las que se dan en la parte inferior.
tt=0:pi/30:2*pi
rho=1:0.1:2
[R,T]=meshgrid(rho,tt)
X=R.*cos(T);
Y=R.*sin(T);
p=10-((R.^2-1).^2./R.^4.*(cos(T)).^2+1./R.^2.*(sqrt(2)-sin(T).*((R.^2+1)./R.^2)).^2);
contour(X,Y,p,500)
axis equal
xlabel('EjeX')
ylabel('EjeY')
