Visuaización de Campos escalares y vectoriales en fluidos (grupo A1)

1 . Campos en un fluido con Obstáculo Circular

Partimos de una superficie plana circular que expresaremos en coordenadas cilíndricas como área de trabajo. O sea, el disco de radio 5 centrado en el origen de coordenadas. Tenemos como obstáculo un círculo de radio 1 centrado en el Origen alrededor del cual observaremos el comportamiento del fluido.

Idea de la situación del fluido.

Malla en polares.

%vectores de las coordenadas
rho=linspace(1,5,50);
phi=linspace(0,2*pi,50);


%malla
h=0.1; %paso
[mrho,mphi]=meshgrid(rho,phi);
figure(1)
X=mrho.*cos(mphi);
Y=mrho.*sin(mphi);
mesh(X,Y,0*X)
axis([-4,4,-4,4])
view(2)

2 . Función Potencial [math]φ[/math] y su Gradiente u

Conocemos una función [math]φ[/math] de la cual sabemos que su gradiente u representa la velocidad de cada punto del fluido.

• [math] φ(ρ,θ,z)=(ρ+ 1/ρ)cosθ [/math]
• Representación de las curvas de nivel de [math]φ[/math] en matlab:

Curvas de nivel de [math]φ[/math].

f = inline('cos(th).*(ro+1./ro)','ro','th');
fro = inline('cos(th).*(1-1./ro.^2)','ro','th');
fth = inline('-(ro.^2+1).*sin(th)./ro.^2','ro','th');


ro = linspace(1,5,10);
th = linspace(0,2*pi,10);


%mallado para las curvas de nivel
[Mro,Mth] = meshgrid(ro,th);
Z=f(Mro,Mth);
X=Mro.*cos(Mth);
Y=Mro.*sin(Mth);

%curvas de nivel
contour(X,Y,Z,12);

• [math]\vec u =∇φ=(u_1,u_2,u_3)= \frac{\partial φ}{\partial ρ}·g^ρ + \frac{\partial φ}{\partial θ}·g^θ + \frac{\partial φ}{\partial z}·g^z [/math]

En coordenadas contravariantes [math] {u_i} [/math] de [math] (g^i) [/math]. Cambiadas a contravariantes [math]{u^i}[/math] de la recíproca [math] (g_i) [/math]:

[math] (u_i)·g^i = (u_i)·G·g_i = \frac{\partial φ}{\partial ρ}·g^1·g^1·g_ρ + \frac{\partial φ}{\partial θ}·g^2·g^2·g_θ + \frac{\partial φ}{\partial z}·g^3·g^3·g_z = (u^1,u^2,u^3)·g_i [/math]

Conociendo los elementos [math]g^i·g^j[/math] de la matriz de gram de la transformación a cilíndricas:

• [math]g^1·g^1=1 [/math]
• [math]g^2·g^2=1/p^2[/math]
• [math]g^3·g^3=1[/math]

entonces:

[math] (u^i)·g_i = [(1-\frac{ 1}{ ρ^2})cosθ]·g_ρ - [(\frac{(p^2+1)}{ p^3})·senθ]·g_θ + 0·g_z [/math]

El programa matlab del campo de velocidades del fluido es:

Campo de velocidades del fluido.

f=inline('cos(th).*(ro+1./ro)','ro','th');

ro = linspace(1,5,80);
th = linspace(0,2*pi,80);

[Mro,Mth] = meshgrid(ro,th);
X = Mro.*cos(Mth);
Y = Mro.*sin(Mth);
fX = (((1./Mro-1./Mro.^3).*cos(Mth)).*X)+(((1./Mro+1./Mro.^3).*sin(Mth)).*Y);
fY = (((1./Mro-1./Mro.^3).*cos(Mth)).*Y)-(((1./Mro+1./Mro.^3).*sin(Mth)).*X);

quiver(X,Y,fX,fY)
axis([-4,4,-4,4]);
view(2)

Si quisiésemos ver que el gradiente es perpendicular a las curvas de nivel en cada punto, realizaremos lo siguiente:

Velocidad ([math]\vec u[/math]) = Gradiente perpendicular a lineas de nive.l

f=inline('cos(th).*(ro+1./ro)','ro','th');

ro = linspace(1,5,20);
th = linspace(0,2*pi,20);

[Mro,Mth] = meshgrid(ro,th);
X = Mro.*cos(Mth);
Y = Mro.*sin(Mth);
fX = (((1./Mro-1./Mro.^3).*cos(Mth)).*X)+(((1./Mro+1./Mro.^3).*sin(Mth)).*Y);
fY = (((1./Mro-1./Mro.^3).*cos(Mth)).*Y)-(((1./Mro+1./Mro.^3).*sin(Mth)).*X);

quiver(X,Y,fX,fY)
axis([-4,4,-4,4]);
view(2)
hold on

%curvas de nivel
Z=f(Mro,Mth);
contour(X,Y,Z,12);

2.1 . Rotacional de u nulo

Sabiendo que se trata de un fluido incompresible podemos observar que el rotacional de [math]u⃗[/math] es nulo:

[math] \nabla\times \vec u⃗ =\frac{ 1}{ ρ}\left| \begin{matrix} \ g_ρ & \ g_θ & \ g_z \\ & & \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ & & \\ u_1 & u_2 & u_3 \end{matrix}\right| [/math]

[math] ∇\times\vec u=(\frac{\partial u^3}{\partial ρ}-\frac{\partial u^2}{\partial θ})·g_ρ-(\frac{\partial u^3}{\partial ρ}-\frac{\partial u^1}{\partial z})·g_θ+(\frac{\partial u^2}{\partial ρ} -\frac{\partial u^1}{\partial θ})·g_z [/math] Las coordenadas en [math]g_ρ[/math] y [math]g_θ[/math] s anulan al derivar, y en la [math]g_z[/math]:

[math] [ -(1-\frac{ 1}{ ρ^2})*senθ + (1-\frac{ 1}{ ρ^2})*senθ ]·g_z [/math] [math] ∇φ=0·g_i [/math]

2.2 . Divergencia de u nula

También vemos que la divergencia se anula:

[math] ∇·\vec u=\frac{ 1}{ √g}·\frac{\partial }{\partial x^i}·(√g·u^i) = [/math]

Sustituyendo:

[math] =\frac{ 1}{ ρ}·[(\frac{\partial }{\partial ρ}·(ρ-\frac{ 1}{ ρ})·cosθ) - (\frac{\partial }{\partial θ}·(1+\frac{ 1}{ ρ^2})·senθ) + 0] = 0 [/math]

2.3 . Módulo de u

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

Fijando [math] ρ = 1[/math] estudiaremos el módulo de la velocidad en la frontera con el obstáculo.

[math] \frac{\partial |\vec u|^2}{\partial θ}=0 [/math] Para calcular máximos y mínimos de la función modulo de [math]|\vec u|^2 [/math]

[math] \frac{\partial |\vec u|^2}{\partial θ}= 8·senθ·cosθ = 0 [/math]

La función del módulo de la velocidad [math]|\vec u|^2 [/math] alcanza puntos críticos en [math] θ = [0, \frac{ \pi}{ 2}, \pi, \frac{ 3·\pi}{ 2}] [/math]

• [math]|\vec u| (0)= 0 [/math] mínimo relativo
• [math]|\vec u| (\frac{ \pi}{ 2})= 2 [/math] máximo relativo
• [math]|\vec u| (\pi)= 0 [/math] mínimo relativo
• [math]|\vec u| (\frac{ 3\pi}{ 2})= 2[/math] máximo relativo

Los puntos de [math] θ = [0,\pi] [/math] tienen velocidad nula y son puntos de remanso.

3 . Dibujar el fluido, lineas de corriente

Existen unas lineas que son tangentes a la velocidad dada por el campo u en cada punto y que determinan la trayectoria de las partículas del fluido. Estas vendrán dadas por el campo v obtenido como:

[math] \vec v=\vec k\times \vec u =ρ\left| \begin{matrix} \ g^ρ & \ g^θ & \ g^z \\ & & \\ 0 & 0 & 1 \\ & & \\ u^1 & u^2 & u^3 \end{matrix}\right| [/math]

[math] \vec k\times \vec u = (-u^2)·g^ρ - (-u^1)·g^θ = (-ρ·u^2,ρ·u^1,0) = v_i·g^i [/math]

• Podemos comprobar que v es irrotacional también:

[math] rot \vec v =1/ρ\left| \begin{matrix} \ g_ρ & \ g_θ & \ g_z \\ & & \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ & & \\ v_1 & v_2 & v_3 \end{matrix}\right| [/math]

donde [math]g_ρ[/math] y [math]g_θ[/math] se anulan al derivar como antes, y [math] [(1+1/ρ^2)cosθ - (1+1/ρ^2)cosθ]· g_z = 0· g_z[/math]

[math]\vec v[/math] es un campo irrotacional también

4 . Líneas de corriente

• Definimos una funcion potencial [math]ψ[/math] de ese campo [math]\vec v[/math]
  

  Campo de velocidades tangente a las lineas de corriente.
• Teniendo en cuenta la definicion de [math]\vec v[/math]:

[math]\frac{\partial ψ}{\partial ρ} = v_1 = -ρ·u^2 = (1+\frac{ 1}{ ρ^2})·senθ [/math]

[math]ψ =\int v_1·dρ =\int (1+\frac{ 1}{ ρ^2})senθ dρ = ρsenθ - \frac{ 1}{ ρ}·senθ + cte(θ)[/math]

[math]\frac{\partial ψ}{\partial θ} = \frac{\partial }{\partial θ}·[ρsenθ - \frac{ 1}{ ρ}·senθ + cte(θ)] = ρcosθ - \frac{ 1}{ ρ}·cosθ + cte'(θ)[/math]

[math]\frac{\partial ψ}{\partial θ} = v_2[/math] [math]ρcosθ - \frac{ 1}{ ρ}·cosθ + cte'(θ) = ρcosθ - \frac{ 1}{ ρ}·cosθ[/math]

[math]cte'(θ)=0[/math] [math]cte[/math] no es [math] cte(θ) [/math] , [math]cte=k [/math]

[math]ψ =ρsenθ - \frac{ 1}{ ρ}·senθ + k = (ρ -\frac{ 1}{ ρ})senθ[/math]

4.1 . Ecuación de Bernoulli

Para liquidos incompresibles tenemos que la ecuación de Bernoulli es:

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

Presiones del fluido.

en campos conservativos como es el caso del fluido a estudiar, [math]g+z = ψ[/math] y [math] |\vec v| [/math] es el vector [math] |\vec u| [/math]. Además nos dicen que [math]p = cte -\frac{1}{2}·ρ_d·|\vec u|^2 [/math]

entonces:

[math]\frac{|\vec u|^2}{2} + ψ + cte' -\frac{ρ_d·|\vec u|^2}{2}=cte[/math]

Particularizando para [math]ρ_d=2[/math]

[math]\frac{|\vec u|^2}{2} + ψ +cte' -\frac {|\vec u|^2}{2}=cte[/math]

donde queda que [math] ψ =cte [/math] cosa que es cierta, cada linea de corriente tiene un potencial constante, verificandose la ecuación de Bernoulli.

5 . Campo u

Siendo el campo [math]\vec u[/math] un campo vectorial gradiente de la función [math]φ[/math], la circulación de este campo conservativo sobre la frontera cerrada del obstáculo circular es nula.

[math]\Gamma= \oint_{C} V \cdot d\mathbf{s}=0[/math]

aplicando esto:

• El teorema de Kutta-Joukowski expresa la fuerza [math]L[/math] como:

[math]L= ρ· \vec u ·\Gamma.[/math]

donde [math]\rho[/math] es la densidad del fluido, [math]\Gamma[/math] la cirlculación del campo [math]V=\vec u[/math] sobre la curva de frontera del obstáculo

[math]L = ρ· \vec u ·0 = 0[/math]

El fluido no ejerce fuerzas sobre el obstáculo.

Ésto ocurre únicamente porque el campo es gradiente y la integral sobre una curva cerrada de estos campos se anula. Es decir tratamos con un fluido ideal conservativo en el que ocurre esta paradoja llamada de d'Alembert.

6 . Cambiando [math] φ(ρ,θ,z)=(ρ+ 1/ρ)cosθ + θ/4\pi[/math]

Aplicando lo mismo que con la otra :

[math]∇φ=(u^i)·g_i = [(1-\frac{ 1}{ ρ^2})cosθ]·g_ρ + [(\frac{1}{ p}+\frac{1}{ p^3})·senθ-\frac{1}{ p^2·4\pi}]·g_θ + 0·g_z[/math]

[math]∇φ=(u_i)·g^i = [(1-\frac{ 1}{ ρ^2})cosθ]·g^ρ + [ \frac{1}{ 4\pi}-(p+\frac{1}{ p})senθ]·g^θ[/math]

• Irrotacional

[math] \nabla\times \vec u =\frac{ 1}{ ρ}\left| \begin{matrix} \ g_ρ & \ g_θ & \ g_z \\ & & \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ & & \\ u_1 & u_2 & u_3 \end{matrix}\right| [/math]

[math] ∇\times\vec u=(\frac{\partial u^2}{\partial ρ} -\frac{\partial u^1}{\partial θ})·g_z = -(1-\frac{ 1}{\ ρ^2})senθ + (1-\frac{ 1}{\ ρ^2})senθ = 0 [/math]

• Divergencia nula

[math] ∇·\vec u = \frac{ 1}{\ ρ}[\frac{\partial }{\partial ρ}·(ρ-\frac{ 1}{\ ρ})cosθ] + \frac{ 1}{\ ρ}[\frac{\partial }{\partial θ}·\frac{ 1}{ ρ·4\pi}-(1+\frac{ 1}{ ρ^2})senθ] = \frac{ 1}{\ ρ}[(1+\frac{ 1}{\ ρ^2})cosθ - (1+\frac{ 1}{\ ρ^2})cosθ] = 0 [/math]

Dibujando en matlab la nueva función, sus curvas de nivel y gradiente:

Gradiente de la nueva [math]\ phi[/math] y sus curvas de nivel.

f=inline('cos(th).*(ro+1./ro)+th./(4*pi)','ro','th');

ro = linspace(1,5,20);
th = linspace(0,2*pi,20);

[Mro,Mth] = meshgrid(ro,th);
X = Mro.*cos(Mth);
Y = Mro.*sin(Mth);
fX = (((1./Mro-1./Mro.^3).*cos(Mth)).*X)+(((1./Mro+1./Mro.^3).*sin(Mth)).*Y)-Y./(4.*pi.*Mro.^2);
fY = (((1./Mro-1./Mro.^3).*cos(Mth)).*Y)-(((1./Mro+1./Mro.^3).*sin(Mth)).*X)+X./(4.*pi.*Mro.^2);

quiver(X,Y,fX,fY)
axis([-4,4,-4,4]);
view(2)
hold on

%curvas de nivel
Z=f(Mro,Mth);
contour(X,Y,Z,20);

Líneas de corriente de la nueva [math]\ phi[/math].

[math]\psi =(ρ+\frac{ 1}{ ρ})·senθ +\frac{ ln ρ}{ 4π} [/math]

Dibujando las nuevas lineas de corriente en matlab

%gradiente de la funcion potencial
f = inline('sin(th).*(ro-1./ro)+log(ro)./(4.*pi)','ro','th');

ro = linspace(1,5,50);
th = linspace(0,2*pi,50);

%mallado para las curvas de nivel
[Mro,Mth] = meshgrid(ro,th);
Z = f(Mro,Mth);
X = Mro.*cos(Mth);
Y = Mro.*sin(Mth);

%dibujo de las curvas de nivel
contour(X,Y,Z,20);