Fluido alrededor de un obstáculo circular (Grupo 20A)

1 Introducción

Un fluido es un medio continuo tal que en reposo no transmite tensión tangencial. No es capaz soportar su propio peso, y tiende a adoptar las características del recipiente en el que se encuentra.

Los fluidos pueden ser líquidos o gases, se van a estudiar los fluidos incompresibles, que son aquellos cuyo volumen no disminuye al ejercer una fuerza sobre él, este tipo de fluidos son los líquidos.

En este artículo se estudiará el comportamiento de estos fluidos al fluir alrededor de un obstáculo con forma circular. Por conveniencia se trabajará con coordenadas cilíndricas.

2 Mallado

Se empezará dibujando un mallado que represente los puntos interiores de la región ocupada por un fluido, que será el exterior del círculo unidad. Para llevar a cabo esta gráfica, se tomará un mallado del anillo comprendido entre los radios 1 y 5 y centro el origen.

La finalidad de realizar un mallado es dividir el dominio continuo (la región ocupada por el fluido) en subdominios más pequeños llamados elementos finitos o celdas. Estos elementos forman una malla que cubre toda la región y permite realizar cálculos y simulaciones numéricas en cada uno de ellos.

Para ilustrar que el fluido ocupa el exterior de un circulo, se dibujarán los ejes en el intervalo [math][-5,5]×[-5,5][/math].

Mallado del obstáculo circular

rho=linspace(1,5,30);                        %Se definen rho y theta 
th=linspace(0,2*pi,30);

[U,V]=meshgrid(rho,th);                      %Se crea el mallado

hold on
X=U.*cos(V);                                 %Se parametriza la superficie
Y=U.*sin(V);
Z=0.*U;
mesh(X,Y,0*Z);                               %Dibujo de la placa 
plot(1*cos(th),1*sin(th),'k','lineWidth',1); %Se representa el obstáculo
axis([-5,5,-5,5]);                           %Se fijan los ejes 
view(2);
title ('Placa');
xlabel 'EJE X'
ylabel 'EJE Y'
axis equal
hold off

3 Función Potencial del Fluido

Una vez definida la región que se va a estudiar, se analizará la velocidad de las partículas del fluido, que viene dada por su función potencial.

La función potencial establecida es [math] \varphi (\rho ,\theta)=(\rho +\frac{1}{\rho})\cos (\theta )- \theta [/math]. Se representará la función potencial, para luego examinar la velocidad de las partículas y el movimiento del fluido.

Mallado del obstáculo circular

rho=linspace(1,5,30);                             %Se define rho y theta
th=linspace(0,2*pi,30);

[U,V]=meshgrid(rho,th);                           %Se crea el mallado

hold on
X=U.*cos(V);                                      %Se parametriza la superficie
Y=U.*sin(V);

f=@(rho,th)(rho+(1./rho)).*cos(th)-th;  %Se define la función potencial
Z=f(U,V);                                         %Se aplica la función potencial
surf(X,Y,Z);                                      %Se representa la función
plot(1*cos(th),1*sin(th),'k','lineWidth',3);      %Se representa 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 del Fluido

La velocidad de las partículas del fluido viene dada por el gradiente de la función potencial. Esto determina el valor y la dirección de máximo crecimiento de la función potencial.

El gradiente de una función en coordenadas cilíndricas se obtiene a partir de la siguiente fórmula: \(\nabla f(\rho, \theta, z) = \frac{\partial f}{\partial \rho} \vec e_\rho + \frac{1}{\rho} \frac{\partial f}{\partial \theta} \vec e_\theta + \frac{\partial f}{\partial z} \vec e_z\)

Siguiendo esta fórmula el gradiente de la función potencial determinada es [math] \vec u=\nabla \varphi=[(1-\frac{1}{\rho^2})\cos(\theta)]\vec e_\rho-[(1+\frac{1}{\rho^2})\sin(\theta) + \frac{1}{\rho}]\vec e_\theta [/math]

Como la representación se realiza en coordenadas cartesianas, se debe pasar el gradiente de la función potencial a estas. Para ello usamos la matriz de cambio de base de cilíndricas a cartesianas y la multiplicamos por el vector dado en coordenadas cilíndricas.

Matriz de cambio de base:

\begin{bmatrix} \cos(\theta) & -\sin(\theta) & 0 \\ \sin(\theta) & \cos(\theta) & 0 \\ 0 & 0 & 1 \end{bmatrix}

Tras el cambio obtenemos;

[math] \vec u=\nabla \varphi=\left( (1 - \frac{1}{{\rho^2}}) \cdot \cos^2(\theta) + (1 + \frac{1}{{\rho^2}}) \cdot \sin^2(\theta) + \frac{sin(\theta)}{\rho}\right) \vec i + \left( (1 - \frac{1}{{\rho^2}}) \cdot \sin(\theta) \cdot \cos(\theta) - (1 + \frac{1}{{\rho^2}}) \cdot \sin(\theta) \cdot \cos(\theta) - \frac{cos(\theta)}{\rho} \right) \vec j [/math]

A continuación, se representa el campo de velocidades.

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

X=U.*cos(V);                                       %Parametrizamos la superficie
Y=U.*sin(V);

f=@(rho,th)(rho+(1./rho)).*cos(th)-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+(sin(V)./rho);           
Cy=(1-(1./U.^2)).*sin(V).*cos(V)-(1+(1./U.^2)).*sin(V).*cos(V)-(cos(V)./rho);
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

Ahora más de cerca comprobamos que \( u = \nabla\varphi \) es ortogonal a las curvas de nivel de \( \varphi \).

4.1 Interpretación y Valor de [math]\vec u [/math]

Utilizando la gráfica anterior, se puede visualizar como el campo de velocidades [math] \vec u=\nabla \varphi [/math], representado por las flechas, es ortogonal en todo punto a las curvas de nivel de [math] \varphi [/math].

Para comprobar matemáticamente esta observación, se elige el vector [math] \vec n [/math] normal a los puntos del obstáculo, y se verifica que cumple la fórmula [math]\vec{u}\cdot \vec{n}=0 [/math]. Esto demuestra que, efectivamente los vectores son ortogonales entre sí.

Se evalúa el gradiente [math] \vec u [/math] en el punto [math]\rho = 1[/math], y por ello se elimina la componente [math]\vec e_\rho [/math] del gradiente. De esta forma, si multiplicamos el gradiente por [math] \vec n [/math], el vector normal, que en este caso es [math]-\vec e_\rho [/math] se tiene que:

[math]-\vec e_\rho \cdot \left ( -2\sin (\theta) -1 \right )\vec e_\theta = 0 [/math]

Se evalúa [math] \vec u [/math] en la frontera del obstáculo:

[math]\vec u(1,\theta,z)=\nabla \varphi(1,\theta,z)= \left (-2\sin(\theta)-1 \right)\vec e_\theta [/math]

4.2 Campo de velocidades lejos del obstáculo

Para poder estudiar la función estando lejos del obstáculo, se hará la siguiente suposición: al estar lejos del obstáculo circular, [math]\rho [/math] es muy grande, y se puede suponer que [math] \frac{1}{\rho } [/math] es despreciable.

Trabajamos sobre la función potencial establecida: [math] \varphi (\rho ,\theta)=(\rho +\frac{1}{\rho})\cos (\theta ) - \theta [/math].

Al despreciar [math] \frac{1}{\rho } [/math] obtenemos lo siguiete; [math] \varphi* = \rho \cos (\theta) - \theta [/math].

Seguidamente calculamos el gradiente de la función potencial hallada:

[math]\nabla f(\rho, \theta, z) = \frac{\partial f}{\partial \rho} \vec e_\rho + \frac{1}{\rho} \frac{\partial f}{\partial \theta} \vec e_\theta + \frac{\partial f}{\partial z} \vec e_z= \nabla \varphi*= [/math] [math]\cos(\theta) \vec{e}_{\rho} - (\sin(\theta) + \frac{1}{\rho}) \vec{e}_{\theta}[/math]

4.3 Rotacional Nulo y Divergencia Nula

El rotacional y la divergencia del campo vectorial proporcionan información sobre las características físicas del fluido.

El operador rotacional indica la dirección y la velocidad del giro del campo en cada punto. Considera el movimiento del fluido y su tendencia a inducir rotación en cada punto.

A continuación se demuestra que el rotacional es nulo:

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

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

[math]\boxed { \nabla \times \bar { u } =0 } [/math]

Obtenemos que el rotacional es 0, por lo que las partículas del fluido no están girando.

La divergencia en cambio, mide la tasa de expansión y contracción de un fluido en cada punto, calculando la diferencia entre el flujo entrante y saliente de la superficie.

A continuación se demuestra que la divergencia es nula:

[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\}= [/math]

[math] = \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\}= [/math] [math] \frac{1}{\rho} \left\{ (1 + \frac{1}{\rho^2})\cos(\theta) -\left(1 + \frac{1}{\rho^2}\right)\cos(\theta) \right\}= 0 [/math]

[math]\boxed { div(\bar { u } )=0 } [/math]

Obtenemos que al igual que el rotacional, la divergencia es también 0, lo que indica que el volumen del fluido es constante, es decir, es un fluido incompresible.

5 Líneas de Corriente

En primer lugar, se calcularán las lineas de corriente del campo [math]\vec{u}[/math], donde las líneas de corriente son aquellas que son tangentes al campo [math]\vec{u}[/math] en cada punto. Dichas líneas son las cuales establecen la trayectoria que sigue las partículas del fluido.

Para su realización, se calculará el campo que es ortogonal en cada punto a [math]\vec{u}[/math]. Este vector se asignará como [math]\vec{v}[/math] y se calculará con la siguiente fórmula: [math]\vec{v}=\vec{k}\times \vec{u}[/math]

En el apartado 4.3, se puede apreciar que el rotacional es nulo y esto implica [math]\vec{v}[/math] sea irrotacional. A su vez, tiene un potencial escalar [math]\psi [/math] el cual tiene [math]\vec{v}[/math] de gradiente. Este potencial es conocido como la función de corriente de [math]\vec{u}[/math].

Para calcular dicho potencial escalar [math]\psi [/math] se resolverá este sistema de ecuaciones: [math]\nabla\psi =\vec{v}[/math] para finalmente dibujar las líneas [math]\psi =cte[/math].

Empezamos las operaciones: [math]\vec{v} =\vec{k}\times \vec{u} = \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) - \frac{1}{\rho}) & 0 \end{vmatrix} = ((1+\frac{1}{\rho^2})\sin{\theta}+\frac{1}{\rho})\vec{e}_\rho + (1-\frac{1}{\rho^2})\cos{\theta}\vec{e}_\theta[/math]

Resolvemos el sistema del potencial sabiendo que:

[math]\frac{\partial \psi}{\partial \rho} = \vec{v}_\rho = ((1+\frac{1}{\rho^2})\sin{\theta}+\frac{1}{\rho}\vec{e}_\rho[/math]

[math]\frac{1}{\rho}\frac{\partial \psi}{\partial \theta} = \vec{v}_\theta = (1-\frac{1}{\rho^2})\cos{\theta}\vec{e}_\theta[/math]

[math]\frac{\partial \psi}{\partial Z} = \vec{v}_Z = 0[/math]

[math]\psi(\rho,\theta,Z) =\int((1+\frac{1}{\rho^2})\sin{\theta} + \frac{1}{\rho})\partial\rho = \sin{\theta}\int(1+\frac{1}{\rho^2})\partial\rho + \int\frac{1}{\rho}\partial\rho= \sin{\theta}(\rho-\frac{1}{\rho})+ log({\rho}) + \delta(\theta,Z)[/math]

[math]\frac{\partial \psi}{\partial \theta} = \cos{\theta}(\rho+\frac{1}{\rho})+\frac{\partial\delta}{\partial\theta}=\rho\cos{\theta}(1-\frac{1}{\rho^2})[/math]

[math]\frac{\partial \delta(\theta,Z)}{\partial \theta} = 0[/math]

[math]\delta(\theta,Z) = cte[/math]

[math]\psi(\rho,\theta,Z) = \log(\rho)+\sin{\theta}(\rho-\frac{1}{\rho})-cte[/math]

Finalmente, obtenemos la ecuación [math]\psi: \log(\rho)+\sin{\theta}(\rho-\frac{1}{\rho})=cte[/math]

6 Puntos de la Frontera con Velocidades Máximas y Mínimas

Seguidamente, donde el módulo de la velocidad es mayor y menor, se calcularán los puntos de la frontera del obstáculo S. Los puntos de remanso, serán en los que la velocidad es nula. En el borde del obstáculo de la gráfica se señalarán los puntos de remanso.

Debido a que [math] \vec{u} =\nabla\varphi [/math] , se tiene que hallar el módulo de la velocidad, [math] \left \| \vec u \right \|[/math]. Este módulo es un escalar, y se analiza en la región S= [math] \left \{(\rho,θ): \rho=1 \right \}[/math]. Al calcular el módulo en dicha región, se obtiene que: [math] \left \| \vec u \right \|= -2\sin(\theta)-1 [/math]

Para encontrar los puntos de inflexión, es decir, los máximos y mínimos, se deriva la ecuación: [math] \frac {\partial}{\partial\theta}\left ( -2sin\theta - 1 \right )=0 \Rightarrow cos(\theta)=0 \Rightarrow \theta=3π/2 [/math] y [math]\theta=π/2 [/math]

Por tanto, obtenemos los siguientes puntos de inflexión; [math] \theta=3π/2 \text{ y } \theta=π/2[/math].

Ya que el módulo de la velocidad es dado en valor absoluto, su valor mínimo posible es [math] \left \| \vec u \right \|=0 \Rightarrow \left \| \vec u \right \|= -2\sin(\theta)-1=0 [/math].

Con estos sacamos los mínimos: [math] \theta=7π/6 \text{ y } \theta=11π/6[/math]

En conclusión, con ayuda del grafico: el máximo absoluto se da en: [math] \theta=π/2 [/math] y los mínimos, es decir los punto de remanso se dan en: [math] \theta=7π/6 \text{ y } \theta=11π/6[/math]

Como cambia el módulo de la velocidad en el borde del obstáculo al aumentar θ

x=linspace(0,2*pi,1000);
y=-2.*sin(x)-1;
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

Para poder calcular los puntos de mayor y menor presión del fluido, se supondrá que la densidad del fluido es constante [math] d = 2 [/math], tal que se verifica la ecuación de Bernoulli: [math] 1/2 d \left \| \vec u \right \|{^2}=cte [/math]

Se le asignará el valor de 10 a la constante anterior, para calcular la presión del fluido.

Aplicando los valores anteriores, la ecuación quedará de la forma: [math] p=10-\left \| \vec u \right \|{^2} [/math]. Debido a que se está observando la región del obstáculo S, se tiene que [math] \left \| \vec u \right \|{^2} [/math] queda en función del parámetro [math] \theta [/math].

[math] p=10- ( -2sinθ - 1 ){^2} = 9 - 4sin{^2}\theta - 4sin\theta[/math]

Para hallar los máximos y mínimos, se deriva la función anterior, y se iguala a 0.

[math]\frac {\partial}{\partial\theta}\left ( 9 - 4sin{^2}\theta - 4sin\theta\right )= -8sin\theta cos\theta - 4cos\theta = 0 \Rightarrow 2sin\theta cos\theta = -cos\theta \Rightarrow sin\theta = -1/2[/math]

De esta forma, se tienen los puntos de inflexión: [math] 7π/6 [/math] y [math] 11π/6 [/math]. Puesto que [math] cos\theta = 0 [/math], puede haber un máximo o mínimo tanto en [math] 3π/2 [/math] como en [math] π/2 [/math].

Los máximos se encuentran en [math] 7π/6 [/math] y [math] 11π/6 [/math], mientras que el mínimo se da en [math] 3π/2 [/math].

A continuación, se proporciona una gráfica que demuestra la presión en el obstáculo circular.

Representación de como cambia la presión en el borde del obstáculo al aumentar θ

tt=linspace(0,2*pi,100)
p=10-(-2.*sin(tt)-1).^2
plot(tt,p)

8 Ecuación de Navier-Stokes

La ecuación de Navier-Stokes es un conjunto de ecuaciones en derivadas parciales, que describen el movimiento de un fluido en función de su campo de velocidades y presiones.

Partiendo de la ecuación de Bernoulli, se busca comprobar que [math] (\vec{u},p) [/math] satisface la ecuación de Navier-Stokes estacionaria, dada por la fórmula:

[math] d(\vec{u}·\nabla)\vec{u} + \nabla p = \mu \Delta \vec{u} [/math]

Debido a que los fluidos incompresibles ocupan siempre el mismo volumen, la viscosidad será nula, por lo que [math] \mu =0 [/math]. Al igual que anteriormente, se considera que la densidad [math] d = 2 [/math]. [math] 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 [/math]

Dados los campos [math] (\vec{u},p) [/math]:

[math] (\vec{u}·\nabla)\vec{u} = u_{j} \frac{\partial u_{i}}{\partial j} \vec{e}_{i} \Rightarrow (\vec{u}·\nabla)\vec{u} = u_{1} \frac{\partial u_{1}}{\partial \rho} \vec{e}_{\rho} + u_{1} \frac{\partial u_{1}}{\partial \theta} \vec{e}_{\theta}+ u_{2} \frac{\partial u_{2}}{\partial \rho} \vec{e}_{\rho}+ u_{2} \frac{\partial u_{2}}{\partial \theta} \vec{e}_{\theta}[/math]

donde el vector [math] \vec{u} = u_{i}\vec{e}_{i} [/math]

De la misma forma, se calcula el gradiente de la ecuación de Bernoulli establecida en el apartado anterior: [math] p = 9 - 4sin{^2}\theta - 4sin\theta [/math]

[math] \nabla p= \frac{\partial p}{\partial \rho} \vec{e}_{\rho} + \frac{\partial p}{\partial \theta} \vec{e}_{\theta} [/math]

[math] \nabla p = [/math]

10 Paradoja de D'Alembert

En primer lugar, p [math] \vec {n} [/math] es lo que se consideraría como la fuerza ejercida por el fluido en los diferentes puntos de la frontera, al ver esto, resulta que el sumatorio de todas las proyecciones de las fuerzas sobre la dirección [math] \vec {i} [/math] da como respuesta que la resultante es nula. En particular, sobre la dirección horizontal del obstáculo el fluido no realiza ningún empuje.

Como ya se ha calculado en el apartado 7, el valor de [math] \ p = 10 - \left\|u\right \|^2 [/math]

Entonces:

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