Visualización de campos escalares y vectoriales en fluidos (Grupo A9)

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en fluidos. Grupo 9-A
Asignatura	Teoría de Campos
Curso	2021-22
Autores	Jesús Doldán Martínez Marina Lancho Mora William Joel Solorzano Farfan Adrian Sebastian Villarroel Zamudio
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción.

Definimos fluido como medio continuo cuyas fuerzas internas se oponen a la interpretación, pero no a la separación, de la materia y no ofrece resistencia al deslizamiento interno mientras este no se produce. Este mismo es incompresible cuando su densidad no es alterada en el tiempo y tiene la capacidad de oponerse a la compresión de si mismo bajo cualquier condición. En el siguiente trabajo expondremos el comportamiento del flujo de un fluido incompresible alrededor de un obstáculo circular mediante la ayuda de sus campos escalares y vectoriales, funciones que permiten modelizar el fluido en cuestión. El campo escalar nos permite asignar a cada punto de nuestro dominio un valor escalar, mientras que el campo vectorial otorga un vector a cada punto del dominio.

2 Representación de los puntos interiores de la región ocupada por un fluido.

Consideramos que nuestro fluido se encuentra en un recinto circular de radio 5 unidades con respecto al origen y que nuestro obstáculo circular es de radio= 1 unidad, también con respecto al centro de coordenadas. Visualizamos este recinto en los ejes definidos, [-5,5] x [-5,5], con ayuda del programa MATLAB. Para ello hemos realizado un mallado circular en coordenadas polares empleando el lenguaje m, obteniendo así el resultado de la figura 1.

%Número de divisiones
n = 50;
%Variables polares
rho = linspace(1,5,n);
tht = linspace(0,2*pi,n);
%Mallado de las polares
[Mrho, Mtht] = meshgrid(rho,tht);
%Mallado pasado a cartesianas
x = Mrho.*cos(Mtht);
y = Mrho.*sin(Mtht);

mesh(x,y,0*x);
axis([-5,5,-5,5]);
view(2);

Figura 1:Malla

3 Función potencial φ

3.1 Campo de velocidades

Figura 2: Campo de velocidad del fluido

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

[math] φ(ρ,θ) = (ρ + \frac{1}{ρ}) \sin(θ)[/math]

Definimos al gradiente en coordenadas polares como:

[math] ∇φ (ρ,θ,z) = (\frac{∂φ}{∂ρ})\vec e_ρ + \frac{1}{ρ}(\frac{∂φ}{∂θ})\vec e_θ + (\frac{∂φ}{∂z})\vec e_z [/math]

Siendo [math](\frac{∂φ}{∂z})\vec e_z = 0[/math] , puesto que la función potencial no depende de la variable [math]z[/math].

El gradiente resultante será:

[math] ∇φ(ρ,θ) = \frac{(ρ^2-1)sin(θ)}{ρ^2}\vec e_ρ + \frac{(ρ^2+1)cos(θ)}{ρ^2}\vec e_θ [/math]

De esta manera obtenemos el campo de velocidades [math] ∇φ [/math] , el cual podemos observar más detenidamente en la figura 2, hecha con el código MATLAB que vemos a continuación.

f = @(rho,tht) (sin(tht).*(rho + 1./rho));
Grho = @(rho,tht) (sin(tht).*(1./rho - 1./(rho.^3)));
Gtht = @(rho,tht) (cos(tht).*(1./rho + 1./(rho.^3)));

n =50;
rho = linspace(1,50,n);
tht = linspace(0,2*pi,n);
[Mrho,Mtht] = meshgrid(rho,tht);

x = Mrho.*cos(Mtht);
y = Mrho.*sin(Mtht);
z = f(Mrho,Mtht);
frho = Grho(Mrho,Mtht);
ftht = Gtht(Mrho,Mtht);
fx = (frho.*x - ftht.*y);
fy = (frho.*y + ftht.*x);
quiver(x,y,fx,fy);
hold on
contour(x,y,z,n)

3.2 Vector gradiente perpendicular a las curvas de nivel de la función potencial

Figura 3: Vectores del campo de velocidad del fluido perpendiculares

En la figura 3 hemos realizado un zoom de la imagen anterior, observando que el vector gradiente es perpendicular a las curvas de nivel de la función potencial, ya que el gradiente nos indica el vector dirección de máxima pendiente de la función en todo su dominio.

3.3 Comprobación

Si denotamos [math]\vec n [/math] como el vector normal a los puntos de la función, podemos comprobar que el producto escalar del [math]\vec n . ∇φ = 0 [/math]

[math]\vec n = \vec e_z [/math]

[math]\vec u = ∇φ = \frac{(ρ^2-1)sin(θ)}{ρ^2}\vec e_ρ + \frac{(ρ^2+1)cos(θ)}{ρ^2}e_θ [/math]

[math]\vec n . \vec u = (0\vec e_ρ + 0\vec e_θ + 1\vec e_z).(\frac{(ρ^2-1)sin(θ)}{ρ^2}\vec e_ρ + \frac{(ρ^2+1)cos(θ)}{ρ^2}\vec e_θ + 0.\vec e_z) = 0 [/math]

Esto muestra que el campo de velocidades es perpendicular a la normal del flujo.

4 Vector [math]\vec u [/math] en coordenadas cartesianas

Suponiendo que [math] \frac{1}{ρ}[/math] es despreciable hacemos el cambio de coordenadas de cilíndricas a cartesianas, obteniendo:

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

También se puede escribir como:

[math]\vec u = \frac{(ρ - \frac{1}{ρ})sin(θ)}{ρ}\vec e_ρ + \frac{(ρ + \frac{1}{ρ})cos(θ)}{ρ}\vec e_θ [/math]

Sabiendo que [math] \frac{1}{ρ}=0[/math] y simplificando queda:

[math]\vec u = sin(θ)\vec e_ρ + cos(θ)\vec e_θ [/math]

Recordamos que para hacer una cambio de coordenadas de cilíndricas a cartesianas también tenemos que cambiar la base:

[math] (\vec e_ρ,\vec e_θ,\vec e_z) → (\vec i,\vec j, \vec k)[/math]
[math] \vec i = (cos(θ)\vec e_ρ - sin(θ)\vec e_θ) [/math]
[math] \vec j = (sin(θ)\vec e_ρ +cos(θ)\vec e_θ) [/math]
[math] \vec k = \vec e_z [/math]

Como podemos observar, el vector [math]\vec u [/math] es el vector [math]\vec j [/math] en coordenadas cartesianas:

[math] \vec u[/math][math](x_1,x_2,x_3)[/math]=[math]\vec j[/math]

5 Rotacional y Divergencia del vector [math]\vec u[/math]

5.1 Rotacional de [math]\vec u[/math]

Definimos el rotacional de un campo vectorial en coordenadas cilindricas como:

[math]∇×\vec u = \frac{1}{ρ}det\begin{pmatrix} \vec{e_ρ} & ρ\vec{e_θ} & \vec{e_z} \\ \frac{∂}{∂ρ} & \frac{∂}{∂θ} & \frac{∂}{∂z} \\ u_ρ & ρu_θ & u_z \end{pmatrix}[/math]

Reemplazando podemos calcular el rotacional:

[math]∇×\vec u = \frac{1}{ρ}det\begin{pmatrix} \vec{e_ρ} & ρ\vec{e_θ} & \vec{e_z} \\ \frac{∂}{∂ρ} & \frac{∂}{∂θ} & \frac{∂}{∂z} \\ \frac{(ρ^2-1)sin(θ)}{ρ^2} & ρ\frac{(ρ^2+1)cos(θ)}{ρ^2} & 0 \end{pmatrix} = \frac{1}{ρ}det\begin{pmatrix} \vec{e_ρ} & ρ\vec{e_θ} & \vec{e_z} \\ \frac{∂}{∂ρ} & \frac{∂}{∂θ} & \frac{∂}{∂z} \\ \frac{(ρ^2-1)sin(θ)}{ρ^2} & \frac{(ρ^2+1)cos(θ)}{ρ} & 0 \end{pmatrix} [/math]
[math] ∇×\vec u = \frac{1}{ρ}([\frac{∂}{∂θ}(0)-\frac{∂}{∂z}(\frac{(ρ^2+1)cos(θ)}{ρ})]\vec e_ρ+[\frac{∂}{∂ρ}(0)-\frac{∂}{∂z}(\frac{(ρ^2-1)sin(θ)}{ρ^2})]ρ\vec e_θ +[\frac{∂}{∂ρ}(\frac{(ρ^2+1)cos(θ)}{ρ})-\frac{∂}{∂θ}(\frac{(ρ^2-1)sin(θ)}{ρ^2})])\vec e_z[/math]
[math] ∇×\vec u = \frac{1}{ρ}([\frac{∂}{∂ρ}(\frac{(ρ^2+1)cos(θ)}{ρ})-\frac{∂}{∂θ}(\frac{(ρ^2-1)sin(θ)}{ρ^2})])\vec e_z = \frac{1}{ρ}(\frac{ρ^2cos(θ)-cos(θ)}{ρ^2}-\frac{ρ^2cos(θ)-cos(θ)}{ρ^2})\vec e_z =0\vec e_z = \vec 0 [/math]

Concluimos que [math]∇×\vec u= \vec 0[/math]

5.2 Divergencia de [math]\vec u[/math]

Definimos divergencia de un de un campo vectorial como:

[math]∇.\vec u =\frac{1}{ρ}(\frac{∂}{∂ρ}(ρu_ρ)+\frac{∂}{∂θ}(u_θ)+\frac{∂}{∂z}(ρu_z))[/math]

Reemplazando obtenemos:

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

Concluimos que [math]∇.\vec u=0[/math]

6 Líneas de corriente del campo [math]\vec u[/math]

Las líneas de corriente son tangentes a cada uno de los puntos de [math]\vec u[/math], es decir, serán ortogonales a un [math]\vec v[/math] que será igual:

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

Definimos como:

[math] \vec k =\vec e_z [/math]
[math] \vec u = \frac{(ρ^2-1)sin(θ)}{ρ^2}\vec e_ρ + \frac{(ρ^2+1)cos(θ)}{ρ^2}\vec e_θ [/math]

Calculando el producto vectorial de ambos vectores:

[math]\vec v[/math] = [math] \vec k × \vec u = det\begin{pmatrix} \vec e_ρ & \vec e_θ & \vec e_z \\ 0 & 0 & 1 \\ \frac{(ρ^2-1)sin(θ)}{ρ^2} & \frac{(ρ^2+1)cos(θ)}{ρ^2} & 0 \end{pmatrix}[/math]
[math]\vec v[/math] = [math] \vec k × \vec u = -\frac{(ρ^2+1)cos(θ)}{ρ^2}\vec e_ρ + \frac{(ρ^2-1)sin(θ)}{ρ^2}\vec e_θ [/math]

Si el vector [math] \vec u [/math] tiene divergencia nula [math] ( \nabla . \vec u [/math]= [math] 0 )[/math] , entonces, el vector [math] \vec v [/math] será irratocional [math] ( \nabla × \vec v [/math]= [math] \vec 0 ) [/math]

Comprobando que [math] \nabla × \vec v [/math]= [math] \vec 0: [/math]

[math] \nabla × \vec v [/math]= [math]\frac{1}{ρ}det\begin{pmatrix} \vec{e_ρ} & ρ\vec{e_θ} & \vec{e_z} \\ \frac{∂}{∂ρ} & \frac{∂}{∂θ} & \frac{∂}{∂z} \\ \frac{-(ρ^2+1)cos(θ)}{ρ^2} & ρ\frac{(ρ^2-1)sin(θ)}{ρ^2} & 0 \end{pmatrix} = \frac{1}{ρ}(\frac{ρ^2sin(θ)-sin(θ)}{ρ^2}-\frac{ρ^2sin(θ)-sin(θ)}{ρ^2})\vec e_z =0\vec e_z = \vec 0 [/math]

Vamos a calcular la función potencial, para ello definimos que:

[math] \nabla ψ = \vec v [/math]
[math] \nabla ψ = \frac{∂ψ}{∂ρ} + \frac{1}{ρ}\frac{∂ψ}{∂θ} + \frac{∂ψ}{∂z} [/math]

Teniendo en cuenta que [math] \frac{∂ψ}{∂z} = 0 [/math]:

[math] \frac{∂ψ}{∂ρ}=\frac{-(ρ^2+1)cos(θ)}{ρ^2} → ψ = -\int (\frac{-(ρ^2+1)cos(θ)}{ρ^2})∂ρ = -\frac{(ρ^2-1)cos(θ)}{ρ} + h(θ)[/math]
[math]\frac{1}{ρ}\frac{∂ψ}{∂θ} = \frac{1}{ρ}(\frac{∂}{∂θ}(\frac{-(ρ^2+1)cos(θ)}{ρ^2}+h(θ)) = \frac{1}{ρ}(\frac{(ρ^2-1)sin(θ)}{ρ}+h'(θ))[/math]

Ahora vamos a igualar [math] \frac{∂ψ}{∂θ}[/math] del vector [math] \vec v [/math] con el [math] \frac{∂ψ}{∂θ}[/math] calculado:

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

Donde podemos deducir que [math] \frac{h'(θ)}{ρ} = 0 → h'(θ) = 0 [/math]

Entonces la función potencial [math] ψ [/math] será:

[math] ψ(ρ,θ) = -\frac{(ρ^2-1)cos(θ)}{ρ} [/math]

Podemos observar el resultado de estas en nuestro campo con la figura 4. Para ello hemos creado este código matlab.

%función
f = @(rho,tht) (sin(tht).*(rho + 1./rho));
Grho = @(rho,tht) (sin(tht).*(1./rho - 1./(rho.^3)));
Gtht = @(rho,tht) (cos(tht).*(1./rho + 1./(rho.^3)));
fc = @(rho,tht) ((-rho + 1./rho).*cos(tht));
%Mallado
n =30;
rho = linspace(1,5,n);
tht = linspace(0,2*pi,n);
[Mrho,Mtht] = meshgrid(rho,tht);
%Coordenadas cartesianas
x = Mrho.*cos(Mtht);
y = Mrho.*sin(Mtht);
z = f(Mrho,Mtht);
z2 = fc(Mrho,Mtht);

frho = Grho(Mrho,Mtht);
ftht = Gtht(Mrho,Mtht);
fx = (frho.*x - ftht.*y);
fy = (frho.*y + ftht.*x);

quiver(x,y,fx,fy);
hold on
contour(x,y,z2,n)

Figura 4: Líneas de corriente del campo [math]\vec u[/math] tangentes

7 Puntos de frontera

Vamos a calcular los puntos de la frontera del obstáculo [math]S[/math] donde el módulo de la velocidad es mayor y menor; y los puntos donde la velocidad es nula (puntos de remanso):

Dando [math] ρ=1 [/math] puesto que obstáculo rodeado por el fluido tiene radio igual a la unidad. Reemplazando en el vector [math] \vec u [/math]:

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

Esta será la dirección donde se encontrarán y valor del campo en todos sus puntos. Podemos concluir que tenemos que saber la variación del ángulo para saber los datos. Por lo que:

[math] -1≤cos(θ)≤1 [/math]

Figura 5: Puntos máximos y mínimos de velocidad del fluido

Tendrá valores máximos cuando:

[math] cos(θ)= -1 → θ = 0 [/math]
[math] cos(θ)=1 → θ = π [/math]

Tendrá valores mínimos cuando:

[math] cos(θ)= 0 → θ ∈ (\frac{π}{2} , \frac{3π}{2}) [/math]

De esta forma, los valores máximos de velocidad serán: [math](1,0)[/math] y [math](-1,0)[/math] y los valores mínimos: [math](0,1)[/math] y [math](0,-1)[/math] como podemos ver en la figura 5 realizada con el siguiente código de lenguaje m.

f = @(rho,tht) (sin(tht).*(rho + 1./rho));
Grho = @(rho,tht) (sin(tht).*(1./rho - 1./(rho.^3)));
Gtht = @(rho,tht) (cos(tht).*(1./rho + 1./(rho.^3)));

n =30;
rho = linspace(1,5,n);
tht = linspace(0,2*pi,n);
[Mrho,Mtht] = meshgrid(rho,tht);

x = Mrho.*cos(Mtht);
y = Mrho.*sin(Mtht);

frho = Grho(Mrho,Mtht);
ftht = Gtht(Mrho,Mtht);
fx = (frho.*x - ftht.*y);
fy = (frho.*y + ftht.*x);

z = sqrt((fx.^2 + fy.^2));

quiver(x,y,fx,fy)
colorbar;
view([45 30]);

Para poder apreciar más la velocidad del fluido de acuerdo a los puntos hemos realizado el siguiente gif, que ha sido conseguido con el siguiente código MATLAB. El gif nos muestra el movimiento del fluido, viendo como evita el obstáculo circular, y empleando los colores para denotar la velocidad.

Grho = @(rho,tht) (sin(tht).*(1 - 1./(rho.^2)));
Gtht = @(rho,tht) (cos(tht).*(1 + 1./(rho.^2)));
fc = @(rho,tht) ((-rho + 1./rho).*cos(tht));

n1 =30;
rho = linspace(1,5,n1);
tht = linspace(0,2*pi,n1);
[Mrho,Mtht] = meshgrid(rho,tht);

x = Mrho.*cos(Mtht);
y = Mrho.*sin(Mtht);
z = fc(Mrho,Mtht);

n = 1000; 
ro_val = 5*ones(n);
tt_val = linspace(pi,2*pi,n);
t = 10;
k = 20;
x_list = zeros(n,t*k);
y_list = zeros(n,t*k);
v_list = zeros(n,t*k);
p_list = zeros(n,t*k);

for i = 1:n
    ro = ro_val(i);
    tt = tt_val(i);
    xl = zeros(1,k*t);
    yl = zeros(1,k*t);
    vl = zeros(1,k*t);
    pl = zeros(1,k*t);
    for j = 1:k*t 
        xl(j) = ro*cos(tt);
        yl(j) = ro*sin(tt);
        vl(j) = sqrt(Grho(ro,tt)^2 + Gtht(ro,tt)^2);
        pl(j) = 10 - vl(j)^2;
        ro = ro + ((Grho(ro,tt))/t);
        tt = tt + 2*(asin(Gtht(ro,tt)/(2*ro)))/t;
    end
    x_list(i,:) = xl;
    y_list(i,:) = yl;
    v_list(i,:) = vl;
    p_list(i,:) = pl;
end    
%Normalizar velocidades
v_max = max(max(v_list));
v_min = min(min(v_list));
v_range = v_max - v_min;
v_norm = (v_list - v_min)./v_range;
%Normalizar presiones
p_max = max(max(p_list));
p_min = min(min(p_list));
p_range = p_max - p_min;
p_norm = (p_list - p_min)./p_range;
%Colores de velocidades
max_c = [0 1 0];
min_c = [0 0 1];
var_c = max_c - min_c;
c_map_v = repmat(var_c,size(v_norm)).*repelem(v_norm,1,3) + repmat(min_c,size(v_norm));
c_map_p = repmat(var_c,size(p_norm)).*repelem(p_norm,1,3) + repmat(min_c,size(p_norm));

sz = 7;
l_range = 1;
for l = 1:k*t
    %contour(x,y,z,n1)
    hold on 
    scatter(x_list(:,l),y_list(:,l),sz,c_map_p(:,l_range:3*l),'o','filled')
    axis([-5 5 -5 5])
    drawnow
    M(l) = getframe;
    l_range = l*3 + 1;
end

Figura 6: Velocidad y comportamiento del fluido

8 Presión del fluido. Ecuación de Bernouilli

La densidad del fluido es constante igual [math] d=2 [/math] y verifica la siguiente ecuación de Bernouilli:

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

Se pide calular la presion siendo la [math] cte=10 [/math]:

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

Recoordar que:

[math] sin^2(θ)+cos^2(θ) = 1 [/math]
[math] cos^2(θ)-sin^2(θ) = cos(2θ) [/math]

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

Reemplezando:

[math] \frac{1}{2}d|\vec u|^2 + p = cte [/math]
[math] \frac{1}{2}(2)(1 + \frac{2cos(2θ)}{ρ^2} + \frac {1}{ρ^4}) + p = 10 → p = 9 - \frac{2cos(2θ)}{ρ^2} - \frac {1}{ρ^4}[/math]

Figura 7: Presión del fluido

Como podemos observar en la figura 6, las presiones máximas y mínimas se alcanzan en los puntos calculados, siendo estos [math](0,1)[/math] y [math](0,-1)[/math], los cuales son opuestos a los puntos de mayor y menor velocidad que veíamos en el apartado anterior acerca de los puntos frontera. Comparando la figura 5 con la figura 7 observamos como son completamente opuestos los valores máximos y mínimos. El código empleado para la realización del último gráfico ha sido el siguiente:

f = @(rho,tht) (sin(tht).*(rho + 1./rho));
Grho = @(rho,tht) (sin(tht).*(1./rho - 1./(rho.^3)));
Gtht = @(rho,tht) (cos(tht).*(1./rho + 1./(rho.^3)));

n =30;
rho = linspace(1,5,n);
tht = linspace(0,2*pi,n);
[Mrho,Mtht] = meshgrid(rho,tht);

x = Mrho.*cos(Mtht);
y = Mrho.*sin(Mtht);

frho = Grho(Mrho,Mtht);
ftht = Gtht(Mrho,Mtht);
fx = (frho.*x - ftht.*y);
fy = (frho.*y + ftht.*x);

z = 10-(fx.^2 + fy.^2);

%forma de encontrar los puntos más altos 
lsz(1) = max(max(z));
lsz(2) = max(max(z-z.*(z == lsz(1))));
k2 =(z >= lsz(2));

%Separa los 2 puntos más altos
xm = reshape(nonzeros(x.*k2),1,[]);
ym = reshape(nonzeros(y.*k2),1,[]);
zm = reshape(nonzeros(z.*k2),1,[]);

surface(x,y,z);
hold on 
scatter3(xm,ym,zm,'ro','filled') %2 puntos altos
colorbar;
view([45 30]);
figure
contour(x,y,z,n)
colorbar;

9 Ecuación de Navier-Stokes estacionaria

Estas ecuaciones pretenden modelar la evolución de las cantidades de la densidad, velocidad y presiones de un fluido, a partir de la Segunda Ley de Newton y la Ley de Conservación de masa. Esta ecuación describe la mecánica de un fluido en general.

Partiendo de la ecuación de Bernouilli:

[math] d(\vec u . \nabla)\vec u + \nabla p[/math] = [math]µ\nabla\vec u [/math]

Vamos a considerar [math] µ = 0 [/math], es decir, viscosidad nula y [math] \nabla . \vec u = 0 [/math] que es la condición de incompresibilidad (volumen constante). Además:

[math] (\vec u . \nabla)\vec u = u_j\frac{\partial u_i}{\partial x_i}\vec e_i[/math]
[math]\nabla\vec u = \nabla u_i . \vec e_i [/math]
[math]\vec u = u_i\vec e_i [/math]

9.1 Demostración de coincidencia

A la ecuación de Navier-Stokes se le aplican las condiciones, queda tal que:

[math] d(\nabla\vec u)\vec u + \nabla\vec u = \vec 0 [/math]

Ecuación de Bernouilli, punto de partida:

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

A la ecuación de Bernuilli aplicamos el operador gradiente a ambos lados

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

Se entiende que [math] \nabla p [/math] se trata del gradiente, que cuya función potencial, es la función p (calculada en el apartado anterior). El gradiente de una constante es [math] \vec 0 [/math].

[math] \frac{1}{2}d\nabla(\vec u\vec u) + \nabla p = \vec 0 [/math]
[math] \frac{1}{2}d(\nabla\vec u\vec u + \vec u\nabla\vec u) + \nabla p= \vec 0 [/math]
[math] \frac{1}{2}d2(\nabla\vec u)\vec u + \nabla p= \vec 0 [/math]
[math] d(\nabla\vec u)\vec u + \nabla p = \vec 0 [/math]

Claramente de puede comprobar que son idénticas, esto se debe principalmente a la similitud de conceptos que emplearon para llegar a su respectivas conclusiones.

10 Variación de presión y velocidad

10.1 Presión

Siendo p:

[math] p = 9 - \frac{2cos(2θ)}{ρ^2} - \frac{1}{ρ^4} [/math]

Dado que la presión variará al cambiar de trayectoria, esta responderá a la forma matemática del gradiente de la presión, cuya fórmula es:

[math] \nabla p = \vec v = (\frac{4cos(2θ)}{ρ^3} + \frac{4}{ρ^5})\vec e_ρ + (\frac{4sen(2θ)}{ρ^3})\vec e_θ [/math]

10.2 Velocidad

Si por ejemplo, la presión sufre un aumento, dado que esta responde o participa en la ecuación de Bernuilli de los apartado anteriores y siendo la igualdad igual a una constante, al aumentar la presión, la velocidad disminuirá tal que:

[math] |\vec u| = \sqrt\frac{2(cte - p)}{d} [/math]

10.3 Corolario

En este caso, la presión será máxima para ρ=1 y para θ=[math] \frac{\pi}{2} y \frac{3\pi}{2} [/math] y por ende la velocidad será mínima en esos puntos. Gráficamente, ese punto corresponde a la zona donde impacta el fluido contra la esfera, en dicha zona se producen elevadas presiones y velocidades mínimas. En las zonas de baja presión, se produce una elevación efectiva del flujo.

11 Teorema de Kutta-Joukowski

Se trata de un teorema básico en la aerodinámica, fue enunciado por Martin Wilhelm Kutta (matemático alemán) y Nikolai Joukowski (ingeniero mecánico ruso) a principios del siglo XX.

El fin del teorema es cuantificar la elevación que se produce por un flujo (en el experimento aire) sobre un cilindro rotando.

Sin embargo, esto produce uan fuerza por unidad de longitud que actúa sobre el cilindro tiene como modulo dG[math] \vec V [/math] y dirección ortogonal a [math] \vec V [/math]

La ecuación matemática es:

[math] Elevación por unidad de volumen = dGV [/math]

Donde d es la densidad del fluido, G es un valor que cuantifica la intensidad del vórtice que se forma (circulación) y V es la velocidad del fluido. Para calcular la intensidad del vórtice se utiliza la siguiente expresión matemática, donde interviene la densidad del fluido de nuevo, ω se trata de la velocidad angular con la que rota en cilindro y r el radio de este.

[math] G = 2dωr^2 [/math]

11.1 Circulación

Se puede llegar a la conclusión de que el flujo/circulación, se puede definir como la integral del campo gradiente [math] \vec u [/math] alrededor de una curva cerrada.

En nuestro caso, se trata de la integral del campo de velocidades del fluido, alrededor de la circunferencias de radio la unidad.

Aplicando el teorema de Stokes, que nos dice que la integral de un campo vectorial continuo y diferenciable a lo largo de una curva cerrada, es igual al flujo del rotacional de una superficie cualquiera de frontera la curva.

[math] \int_{L} \vec u \cdot d(\vec r) = \int_{S} (\nabla \times \vec u) \vec n \cdot dS [/math]

Tal y como se ga calculado en el apartado 5.1 el rotacional de este campo vectorial es nulo, por lo que la circulación será nula.

[math] \int_{S} (\nabla \times \vec u) \vec n \cdot dS = \int_{S} (0)\vec n \cdot dS = 0 [/math]

Además, el rotacional, del gradiente de un función potencial es nulo.

[math] \nabla \times (\nabla ∂) = 0 [/math]

11.2 Paradoja de D'Alambert

La paradoja de Jean Le Rond D'Alambert (matemático francés) es una contradicción matemática a la que él llego. Además, fue de los primeros en emplear la formulación del flujo potencial. Esta paradoja de basa en que no actuarán fuerzas sobre un cuerpo que se mueve a velocidad constante a través de una masa de fluido incompresible no viscoso. Esto contradice a la observación, de ahí la contradicción.

12 Curvas de presión

A continuación, se adjunta el algoritmo en lenguaje M y las respectivas curvas nivel que este proporciona.

%La función y sus gradientes en polares
f = @(rho,tht) (sin(tht).*(rho + 1./rho));
Grho = @(rho,tht) (sin(tht).*(1 - 1./(rho.^2)));
Gtht = @(rho,tht) (cos(tht).*(1 + 1./(rho.^2)));
%Mallado de la región
n =30;
rho = linspace(1,10,n);
tht = linspace(0,2*pi,n);
[Mrho,Mtht] = meshgrid(rho,tht);
%A cartesianas
x = Mrho.*cos(Mtht);
y = Mrho.*sin(Mtht);
%Gradiente en polares
frho = Grho(Mrho,Mtht);
ftht = Gtht(Mrho,Mtht);
%Cambio del gradiente a cartesianas
fx = (frho.*x - ftht.*y)./sqrt(x.^2 + y.^2);
fy = (frho.*y + ftht.*x)./sqrt(x.^2 + y.^2);
%Velocidad en el punto dado por vectores
v = sqrt(fx.^2 + fy.^2);
z = 10-v.^2;

%Gráfica de las curvas de presiones
contour(x,y,z,n)
axis([-5 5 -5 5])

Figura 8: Curvas de presión del fluido

Tal y como se puede comprobar, en las zonas en las que donde la velocidad es menor, hay un considerable aumento de la presión.

Se puede comprobar como las líneas de flujo más próximas al cilindro van reduciendo su drásticamente su velocidad conforme se acercan al cilindro de radio 1, aumentado la presión. Para posteriormente, incrementar la velocidad mientras rodea el obstáculo, alcanzando su punto máximo en ρ=1 y θ=[math] 0 [/math] y [math] \pi [/math] con valores que pueden alcanzar el doble que en puntos más alejados del cilindro, como ρ=4 y θ=[math] 0 [/math].