Flujo de un fluido incompresible alrededor de un obstáculo circular (Grupo A)

1 Introducción

Un fluido se define como una sustancia en la que sus moléculas sufren una atracción débil entre si, lo que se llaman fuerzas cohesivas. Un fluido incompresible es aquel en el que las fuerzas cohesivas son suficiente fuertes para mantener el volumen constante,pero lo suficientemente débiles como para permitir una total deformabilidad.

En este informe se estudiará el flujo de un fluido incompresible alrededor de un obstáculo. Para ello se utilizarán dos funciones que modelizarán el comportamiento del fluido de la manera que se describe a continuación.

En primer lugar parametrizamos con coordenadas cilíndricas el anillo comprendido entre las circunferencias de radios 1 y 5. A partir de esta parametrización se obtendrá el mallado:

Mallado de un fluido incompresible alrededor del anillo.

Este mallado es el que se usará para el resto de programas que se utilizan en este informe. Se informa que algunos de los programas sucesivos no incluyen la creación de este mallado y por lo tanto, para su funcionamiento necesitan la ejecución previa de el programa que crea el mallado descrito. El programa para crear y representar el mallado es el siguiente:

% Mallado
h=0.1;
th= 0:h:2*pi+h;
p= 1:h:5;
[P,Th]= meshgrid (p,th);
X= P.*cos (Th);
Y= P.*sin (Th);
mesh(X,Y,0.*X);
axis([-4,4,-4,4]);
view(2);

2 Velocidad de las partículas del fluido

La velocidad de las partículas de un fluido viene dada por el gradiente de la función potencial:

                 [math]\varphi(\rho,\theta)=(\rho+1/\rho)\cos \theta [/math]

Curvas de nivel de la función potencial

La representación de la función potencial se obtiene a través del programa:

%Representacion función potencial
Z=(P+1./P).*cos(Th);
plot3(X,Y,Z);
axis image
xlabel('Eje x')
ylabel('Eje y')
zlabel('Eje z')

Tras esto calculamos el campo ([math]\vec u[/math]) que representa las velocidades de las partículas del fluido. Para ello debemos hacer el gradiente de un campo escalar. La fórmula general es la siguiente: [math]∇f=\frac{df}{dx^i }(\vec g^i )[/math]

Por lo tanto: [math]\vec u=\nabla \varphi= cosθ(1 - \frac{1}{ρ^2})\vec g^ρ -sin⁡θ (ρ+ \frac{1}{ρ})\vec g^θ [/math] Esto está en coordenadas polares, para representarlo en Octave debemos a pasarlo a coordenadas paramétricas. [math]\vec u=(cosθ - \frac{⁡cos⁡θ}{ρ^2})(cos⁡θ \vec i +sin⁡θ \vec j)-sin⁡θ ( ρ+\frac{1}{ρ})(-\frac{sin⁡θ}{ρ} \vec i + \frac{cos⁡θ}{ρ} \vec j)[/math]

El resultado de representar el campo vectorial:

Campo vectorial sobre las curvas de nivel de la función potencial

Y el programa hecho para ello es:

%% Curvas de nivel
contour(X,Y,Z,100)
hold on

%% Gradiente
fp=cos(Th).*(1-1./P.^2);
fth=-(1./P+1./P.^3).*sin(Th);
U=fp.*cos(Th)-fth.*P.*sin(Th);
V=fp.*sin(Th)+fth.*P.*cos(Th);
quiver(X,Y,U,V)
axis equal
axis image
xlabel('Eje x')
ylabel('Eje y')
hold off 

Se observa que el campo [math]\vec u[/math] es perpendicular a las curvas de nivel de la función potencial.

Perpendicularidad del campo vectorial a las curvas de nivel de la función potencial

Esto es algo que se cumple para todo campo [math]\vec u[/math] con potencial escalar. A continuación se verificará que la divergencia y rotacional del campo [math]\vec u[/math] son nulos. La divergencia de un campo vectorial se define como la traza del tensor [math]\ ∇(\vec u )[/math]. Para un campo expresado en la base natural g: [math]\ ∇\vec u= \frac{1}{√g} \frac{d(√g u^i)}{dx^i}[/math] Por lo tanto reescribimos el campo de la forma [math]\vec u =u^i (\vec g_i )[/math] esto se hará aplicando la matriz G de las coordenadas cilíndricas. [math] \vec g^i=G^{ij} \vec g_i [/math]

Por lo tanto: [math] \vec u=( cos⁡θ- \frac{cos⁡θ}{ρ^2})\vec g_ρ - \frac{1}{ρ^2}(ρ+ \frac{1}{ρ})sin⁡θ \vec g_θ [/math] Y ahora se calcula la divergencia [math] ∇ \vec u=\frac{1}{ρ} (\frac{(∂(ρ(cosθ - \frac{cosθ}{ρ^2 })}{∂ρ} + \frac{∂(\frac{1}{ρ} (-ρsin⁡θ- \frac{sinθ}{ρ})}{∂θ})= cosθ(\frac{1}{ρ}+\frac{1}{ρ^3}) -cos⁡θ (\frac{1}{ρ}+ \frac{1}{ρ^3})=0 [/math]

Con lo que queda demostrado que la divergencia es nula. Ahora se procederá al cálculo del rotacional. Se define rotacional como el doble vector axial asociado a la pare antisimétrica de [math] ∇\vec u [/math]. Si [math]\vec u [/math] viene dado en una base natural cualquiera:

Con lo que queda demostrado que el rotacional es nulo. Que la divergencia del campo sea nula está asociada en este caso al hecho de que localmente el fluido mantiene su volumen, es decir, se trata de la condición de incompresibilidad. Que el rotacional sea nulo significa que el giro del campo respecto a un punto es nulo. Esto último puede parecer contradictorio, pues el fluido ha de girar en torno al obstáculo para sortearlo. Esta aparente contradicción se resolverá más adelante.

3 Líneas de Corriente

Ahora se desea conocer las líneas de corriente del campo [math] \vec u [/math], es decir, las trayectorias que siguen las partículas del fluido. Estas son las líneas tangentes al vector en cada punto. Éstas se obtendrán calculando el potencial del campo ortogonal a [math] \vec u [/math]. El campo ortogonal a [math] \vec u [/math] se define como: [math] \vec v=\vec k × \vec u [/math] Para realizar este cálculo expresaremos el campo en coordenadas covariantes: [math] \vec v_i = G_{ij} \vec v^j [/math] [math] \vec u= (cos⁡θ-\frac{cos⁡θ}{ρ^2})\vec g^ρ -sin⁡θ(ρ+\frac{1}{ρ})\vec g^θ [/math]

[math] \vec k=\vec g^z [/math]

[math] \vec v =\vec k × \vec u= \frac{1}{ρ}(cos⁡θ-\frac{cos⁡θ}{ρ^2})\vec g_θ + sin⁡θ(1+\frac{1}{ρ^2})\vec g_ρ [/math]

[math] \vec v = sin⁡θ(1+\frac{1}{ρ^2})\vec g^ρ + cos⁡θ(ρ-\frac{1}{ρ})\vec g^θ [/math]

Una vez obtenido, comprobamos que el campo es irrotacional, pues debe serlo ya que [math] \vec u [/math] es de divergencia nula:

Efectivamente se comprueba que [math] \vec v [/math] es irrotacional, y como la región en la que está definida es simplemente conexa, el campo tendrá potencial escalar, ψ, conocida como función corriente de [math] \vec u [/math]:

[math] \frac{∂ψ}{∂ρ}=sin⁡θ(1+\frac{1}{ρ^2}) [/math]

[math] \frac{∂ψ}{∂θ}=cos⁡θ(ρ-\frac{1}{ρ})[/math]

[math] \frac{∂ψ}{∂z}=0 [/math]

[math] ψ_ρ=∫sinθ(1+\frac{1}{ρ^2})∂ρ=sin⁡θ(ρ-\frac{1}{ρ})+f(θ,z) [/math]

[math] ψ_θ=∫cosθ(ρ-\frac{1}{ρ})∂θ= sin⁡θ(ρ-\frac{1}{ρ})+g(z)[/math]

[math] ψ_z=0 [/math]

[math]ψ= sin⁡θ(ρ-\frac{1}{ρ})[/math] 

Las líneas de corriente quedan definidas por las curvas de nivel de la función de corriente de [math] \vec u [/math], es decir, ψ=cte

Líneas de corriente en las que el potencial escalar es constante

Para esta representación se ha utilizado un programa similar al usado para representar φ:

%% Mallado
th= 0:0.1:2*pi+0.1;
p= 1:0.1:5;
[P,Th]= meshgrid(p,th);
X= P.*cos(Th);
Y= P.*sin(Th);
Z=sin(Th).*(P-1./P);

%% Líneas de corriente
contour(X,Y,Z,100)

Si además representamos sobre éstas el campo [math] \vec u [/math], se observa que son tangentes.

4 Relación entre velocidad y presión del fluido

A continuación se plantea cuáles serán los puntos de la frontera del obstáculo que tienen mayor y menor velocidad. Para conocer estos puntos habrá que hallar los máximos y mínimos de la función [math]|\vec u |[/math] para ρ=1. El tomar otra constante lo único que hará será desplazar verticalmente la función.

[math] |\vec u|=√((cos⁡θ(1-\frac{1}{ρ}))^2+(-sinθ⁡(ρ+\frac{1}{ρ}))^2\frac{1}{ρ^2} [/math]

Para facilitar la derivación se usará la función[math] |\vec u|^2 [/math] pues lo máximos y mínimos se mantienen.

Para ρ=1; [math] |\vec u|^2=4sin^2⁡θ [/math]

Por lo tanto [math] ∇|\vec u|^2= 8sin⁡θcos⁡θ \vec g_θ [/math]

Igualando a cero se observa que habrá puntos críticos en [math] θ=0,\frac{π}{2},π,\frac{3π}{4} [/math] Observando la representación de [math] \vec u [/math] se observa que [math] \vec u [/math](1,0) [math] \vec u (1,π) [/math] son mínimos y que [math] \vec u (1,\frac{π}{2})[/math] y [math] \vec u (1,\frac{3π}{4})[/math] son máximos.

Módulo de la velocidad

Módulo de la velocidad

Los programas utilizados son:

%% Representación módulo de la velocidad
Vel=sqrt(((sin(Th).*(1+1./P.^2)).^2+(cos(Th).*(1-1./P.^2)).^2));
plot3(X,Y,Vel);
axis image
xlabel('Eje x')
ylabel('Eje y')
zlabel('Eje z')

%% Representación módulo de la velocidad 2
h=pcolor(X,Y,Vel);
set(h,'EdgeColor','none');

Ahora se supondrá que la densidad del fluido es constante, [math] ρ_d=2 [/math] y que se verifica la ecuación de Bernoulli:

[math] \frac{1}{2} ρ_d |\vec v|^2+p=cte [/math]

Se tomará la constante 1. Por lo tanto [math] p=1-|\vec v|^2=1-((sin⁡θ (1+\frac{1}{ρ^2 })^2+(cos⁡θ(1-\frac{1}{ρ^2})^2 ) [/math]

Densidad del fluido para cte=1

Densidad del fluido para cte=1

Para la primera representación se ha utilizado el siguiente programa:

%Representacion presiones
Pr=1-((sin(Th).*(1+1./P.^2)).^2+(cos(Th).*(1-1./P.^2)).^2);
plot3(X,Y,Pr);
axis image
xlabel('Eje x')
ylabel('Eje y')
zlabel('Eje z')

Y para la segunda:
%Represenatción presiones 2
h=pcolor(X,Y,Pr);
set(h,'EdgeColor','none');

Se observa que la presión es máxima donde la velocidad es mínima y la presión es mínima donde la velocidad es máxima, tal y como indica la ecuación de Bernoulli. Observando las gráficas realizadas podemos ver como variará la velocidad y presión de una partícula del fluido al rodear el obstáculo: Al acercarse la partícula al obstáculo su velocidad irá disminuyendo hasta alcanzar su mínimo y la presión aumentando hasta alcanzar su máximo. La partícula que choca justo con el centro del obstáculo llegará a detenerse. Tras alcanzar este mínimo en la velocidad, la partícula se irá acelerando al sortear el obstáculo alcanzando su máxima velocidad al llegar al eje de simetría paralelo al eje Y; hasta este momento, la presión ha ido disminuyendo hasta alcanzar su mínimo, la partícula se irá frenando hasta alcanzar de nuevo el mínimo anteriormente descrito. En consecuencia su presión aumentará hasta alcanzar el máximo anterior. Por último la partícula se acelerará según se aleja del obstáculo y su presión disminuirá hasta alcanzar las condiciones de velocidad y presión previas a la acción del obstáculo.

Tras este análisis podemos resolver la aparente contradicción de la que se hablaba antes. Como se había dicho, viendo las lineas de corriente se observa que el fluido gira en torno al obstáculo para esquivarlo, sin embargo el rotacional del campo de velocidades es nulo. Esto se debe a las diferencias de velocidades existentes en las trayectorias en torno al obstáculo, lo cual "compensa" el giro. En consecuencia, si un objeto fuera arrastrado por este fluido recorrería la trayectoria correspondiente a las lineas de corriente sin girar con respecto a si mismo.

5 Fuerza ejercida sobre el obstáculo

Como [math] \vec u [/math] es conservativo y γ una curva cerrada: [math] ∫ \vec u \vec {dr} = ∫_0^{2π} \vec u \vec {dr}= φ(2π)- φ(0)=0 [/math] Comprobándose así el Teorema de Kutta-Joukowski que relaciona la fuerza de sustentación generada por un cilindro recto con la velocidad del fluido a lo largo del cilindro; El teorema se refiere al flujo bidimensional alrededor de un cilindro (o un cilindro de envergadura de ala infinito) y determina la sustentación generada por unidad de envergadura.

En este caso se relaciona el teorema con la fuerza que ejerce el fluido sobre el obstáculo (círculo unidad), proporcional a esta circulación; que no es más que la integral de línea de la velocidad del fluido, en una curva cerrada que contiene al obstáculo. Al ser [math] \vec u [/math] conservativo y γ una curva cerrada la circulación es nula; por tanto, el fluido no ejerce ninguna fuerza sobre el obstáculo, en contra de la intuición. Esto es lo que se conoce como la paradoja de D'Alembert.

6 Modelo matemático alternativo del flujo

Debido al resultado paradójico obtenido se realizará el análisis de una nueva función para modelizar el flujo de un fluido. Esta función será: [math] φ(ρ,θ)=(ρ+\frac{1}{ρ})cosθ+ \frac{θ}{4π} [/math]

6.1 Función potencial

Vamos a calcular la velocidad del fluido para una función potencial diferente: [math] φ(ρ,θ)=(ρ+\frac{1}{ρ})cosθ+ \frac{θ}{4π} [/math]

Velocidad de la nueva función potencial del fluido

En la representación de su gráfica se observa que es muy similar a la hallada anteriormente pero, debido al sumando que depende de θ, la función presenta una discontinuidad. Al ser discontinua en [math] R^2 [/math], no puede ser escalar potencial en esa región. Por lo tanto debemos definir una región en la que si sea continua. Esta región será [math] R^2 [/math] excluyendo θ=0.

La representación ha sido obtenida a través del siguiente programa:

% Mallado
h=0.1;
th= 0:h:2*pi+h;
p= 1:h:5;
[P,Th]= meshgrid (p,th);
X= P.*cos (Th);
Y= P.*sin (Th);
mesh(X,Y,0.*X);
axis([-4,4,-4,4]);
view(2);

%Representacion función potencial
Z=(P+1./P).*cos(Th)+Th./(4.*pi);
plot3(X,Y,Z);
axis image
xlabel('Eje x')
ylabel('Eje y')
zlabel('Eje z')

6.2 Velocidad

Hallamos el campo de velocidades del nuevo fluido a través del gradiente de la función anterior: [math] \vec u = ∇φ=cosθ(1 - \frac{1}{ρ^2}) \vec g^ρ+(-sin⁡θ(ρ+\frac{1}{ρ}))+\frac{1}{4π}\vec g^θ ) [/math]

La representación del campo vectorial sobre las curvas de nivel de la función potencial se realizará en todo [math] R^2 [/math], teniendo en cuenta que no es conservativo en esta región.

Campo vectorial sobre las curvas de nivel de la función potencial

Igual que en el otro fluido, el campo [math] \vec u [/math] es perpendicular a las curvas de nivel.

%% Curvas de nivel
contour(X,Y,Z,100)
hold on
%% Gradiente
fp=cos(Th).*(1-1./(P.^2));
fth=-(1./P+1./P.^3).*sin(Th)+1./(4.*(P.^2).*pi);
U=fp.*cos(Th)-fth.*P.*sin(Th);
V=fp.*sin(Th)+fth.*P.*cos(Th);
quiver(X,Y,U,V)
axis equal
axis image
xlabel('Eje x')
ylabel('Eje y')
hold off

Estudiamos su divergencia y rotacional del campo [math] \vec u [/math]. En coordenadas contravariantes: [math] \vec u= (cosθ-\frac{cos⁡θ}{ρ^2}\vec ) g_ρ - (\frac{1}{ρ^2})(ρ+\frac{1}{ρ}sinθ ) \vec g_θ [/math]

Ahora comprobamos que la divergencia es nula:

[math] ∇ \vec u= \frac{1}{ρ}(\frac{(∂(ρ(cosθ-\frac{cosθ}{ρ^2 })}{∂ρ})+\frac{∂(-sinθ-\frac{sinθ}{ρ^2} + \frac{1}{4ρπ} )}{∂θ})= cosθ(\frac{1}{ρ}+ \frac{1}{ρ^3})-cosθ(\frac{1}{ρ} + \frac{1}{ρ^3})=0 [/math]

Asímismo comprobamos el rotacional:

Que el rotacional sea nulo es una condición necesaria, pero no suficiente para que un campo sea conservativo. Por ello, este resultado no contradice lo anteriormente dicho, [math] \vec u [/math] no será conservativo en esta región.

6.3 Líneas de corriente

Para conocer la trayectoria de las nuevas líneas de corriente del campo [math] \vec u [/math], calculamos el potencial del campo ortogonal a [math] \vec u [/math] .

                                    [math] \vec v= \vec k × \vec u [/math]

Expresamos el campo en coordenadas covariantes. [math] \vec v =\vec k × \vec u= cosθ(\frac{1}{ρ^1} + \frac{1}{ρ^3})\vec g_θ + (sinθ(1+\frac{1}{ρ^2} )- \frac{1}{4πρ})\vec g_ρ [/math]

[math] \vec v =(-sinθ(1+\frac{1}{ρ^2} )-\frac{1}{4πρ})\vec g^ρ + cosθ(ρ-\frac{1}{ρ}) \vec g^θ [/math]

Una vez obtenido, comprobamos que el campo es irrotacional, pues debe serlo ya que [math] \vec u [/math] es de divergencia nula:

Efectivamente se comprueba que [math] \vec v [/math] es irrotacional, y como la región en la que está definida es simplemente conexa, el campo tendrá potencial escalar, ψ, conocida como función corriente de [math] \vec u [/math]:

[math] \frac{∂ψ}{∂ρ}=sinθ(1+\frac{1}{ρ^2})-\frac{1}{4π} ; \frac{∂ψ}{∂θ}=cosθ(ρ-\frac{1}{ρ}) ; \frac{∂ψ}{∂z}=0 [/math]

                             [math] ψ= sinθ(ρ-\frac{1}{ρ})-\frac{1}{4π}ln ρ [/math]

Las líneas de corriente quedan definidas por la función potencial, ψ=cte

%% Mallado
th= 0:0.1:2*pi+0.1;
p= 1:0.1:5;
[P,Th]= meshgrid(p,th);
X= P.*cos(Th);
Y= P.*sin(Th);
Z=sin(Th).*(P-1./P)-(1/(4*pi))*log(P);

%% Líneas de corriente
contour(X,Y,Z,100)

  Las líneas de corriente quedan definidas por las curvas de nivel de la función de corriente de [math] \vec u [/math], ψ=cte.

6.4 Fuerza ejercida sobre el obstáculo

Como el campo no es conservativo en la región que incluye la curva que representa el obstáculo, la circulación, a través de dicha curva, no será nula. Significando, según el Teorema de Kutta-Jukowski anteriormente explicado, que el fluido ejercerá una fuerza no nula sobre el obstáculo.

Por lo tanto calcularemos la circulación del campo a través de la curva, definida por ρ=1; θ=t; variando t entre 0 y 2π

[math] ∫ \vec u= ∫_0^{2π}(-2sint+\frac{1}{4π})\frac{∂θ}{∂t}dt= ∫_0^{2π}-2sint+\frac{1}{4π}dt= \frac{1}{2} [/math]

Este es un resultado más acorde con la realidad, por lo tanto, podemos decir que este modelo matemático se ajusta mejor al flujo real de un fluido incompresible alrededor de un obstáculo circular.

7 Autores

Grupo A:

• Juan Junguito
• Paula Costa Palmero
• Álvaro Roales Blanco
• Rocío López de Lucas