Visualización de campos escalares y vectoriales en fluidos. (Grupo 23-C)

1 Introducción

Vamos a estudiar el flujo de un fluido incompresible alrededor de un obstáculo con forma circular de radio igual a 1, ósea un círculo unidad. Para la realización de este trabajo usaremos tanto coordendas cartesianas como con coordenadas Cilíndricas (polares).

2 Mallado

La región ocupada por el fluido será el exterior del círculo unidad tal y como indica a continuación el mallado que hemos dibujado. Hemos tomado un mallado del anillo comprendido entre los radio 1 y 5 y centro el origen, con unos ejes dentro del intervalo [-4,4]x[-4,4].

El código matlab con el que lo hemos realizado es el siguiente:

u=linspace(1,5,100);    % Creamos el vector u en el intervalo[1,5]
v=linspace(0,2*pi,100); % Creamos el vector v en el intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v);  % matrices de las coordenadas u y v
xx=uu.*cos(vv);         % Parametrizamos en coordendas cartesianas y así definimos el anillo
yy=uu.*sin(vv);
mesh(xx,yy,0*xx)        % Dibujamos el mallado
hold on
x=cos(v);
y=sin(v);
plot(x,y,'k','Linewidth',1) % Dibujamos el obstáculo: círculo unidad
axis([-4,4,-4,4])           % Definimos los ejes en  [-4,4]x[-4,4]
view(2) % Visualizamos lo dibujado en el plano x e y

Mallado de la región

3 Velocidad de las partículas del fluido

La velocidad de las partículas del fluido viene dada por un campo vectorial [math]\vec u [/math] cuya función potencial es φ=(ρ+1/p)*cos⁡θ. Luego [math]\vec u[/math] es el gradiente de la función potencial: : [math]\vec u=∇φ[/math] Realizaremos los cálculos tanto en cartesianas como en cilíndricas:  : [math]φ(ρ,θ)=(ρ+1/p)*cos⁡θ = φ(x,y)=x+\frac{x}{x^2+y^2}[/math] En el caso de Cartesianas:: [math]φ(x,y)=x+\frac{x}{x^2+y^2}[/math]

[math]\vec u=∇φ(x,y)=\frac{∂φ}{∂x} \vec i +\frac{∂φ}{∂y} \vec j =(1+\frac{(y^2-x^2)}{(x^2+y^2 )^2 }) \vec i -\frac{2xy}{(x^2+y^2 )^2} \vec j[/math]

En el caso de Cilíndricas calculamos tanto las coordenadas contravariantes como las covariantes: : [math]φ(ρ,θ)=(ρ+1/p)*cos⁡θ[/math]

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

Las representaciones de estos dos campos, el campo vectorial [math]\vec u [/math] y el campo escalar φ, función potencial de [math]\vec u [/math], son las siguientes:

%Representación del Campo escalar φ
u=linspace(1,5,100);
v=linspace(0,2*pi,100); 
[uu,vv]=meshgrid(u,v); %mallado sobre el que representar
xx=uu.*cos(vv);       
yy=uu.*sin(vv);
fi=inline('x+x./(x.^2+y.^2)','x','y'); %Definimos la función potencial fi
f=fi(xx,yy);
contour(xx,yy,f,30) %Dibujamos sus líneas de nivel
axis([-4,4,-4,4])      
hold on
x=cos(v);
y=sin(v);
plot(x,y,'k','Linewidth',1) % Dibujamos el campo escalar fi
axis([-4,4,-4,4])  % Seleccionamos la región donde lo dibuja
axis equal
 
%Representación del Campo vectorial u 
u=linspace(1,5,15);
v=linspace(0,2*pi,15); 
[uu,vv]=meshgrid(u,v); %mallado sobre el que representar
xx=uu.*cos(vv);       
yy=uu.*sin(vv);
fx=inline('1+(y.^2-x.^2)./((y.^2+x.^2).^2)','x','y');    %Definimos la función del campo u de velocidades en x 
fy=inline('-2*(x.*y)./((y.^2+x.^2).^2)','x','y');        %Definimos la función del campo u de velocidades en y 
u1=fx(xx,yy);
u2=fy(xx,yy);
quiver(xx,yy,u1,u2)     %dibujamos el campo vectorial u
axis([-4,4,-4,4])
view(2) %Para visualizar lo dibujado

hold off

Campo vectorial [math]\vec u [/math] ortogonal a las curvas de nivel de φ

Como podemos observar, el campo vectorial [math]\vec u [/math] de las velocidades de las partículas es ortogonal a las curvas de nivel de φ ; porque al ser el gradiente de la función potencial representa la línea de máxima pendiente del campo escalar φ.

Campo [math]\vec u [/math] ortogonal a las líneas de nivel de φ zoom

4 Rotacional de [math]\vec u[/math] nulo

El rotacional de un campo nos indica la tendencia de este a inducir rotación alrededor de un punto.

El rotacional de nuestro campo [math]\vec u[/math] es nulo ya que todo campo potencial, es decir, expresable como el gradiente de un potencial escalar es irrotacional.

Vamos a comprobarlo entonces analíticamente realizando los correspondientes cálculo en coordenadas cilíndricas:

[math]\nabla\times\vec{u}= \frac {1}{g^\frac{1}{2}}\begin{pmatrix} \vec{g_ρ} & \vec{g_θ} & \vec{g_z} \\ \frac{ \partial}{\partial ρ} & \frac{\partial }{\partial θ} & \frac{\partial }{\partial z}\\ u_1 & u_2 & u_3 \end{pmatrix}= \frac {1}{ρ} \begin{pmatrix} \vec{g_ρ} & \vec{g_θ} & \vec{g_z} \\ \frac{ \partial}{\partial ρ} & \frac{\partial }{\partial θ} & \frac{\partial }{\partial z}\\ cos(θ)(1 - \frac {1}{ρ^2})& - sen (θ) (ρ + \frac {1}{ρ}) & 0 \end{pmatrix}=[/math] [math]=\frac {1}{ρ}(\frac{\partial}{\partial z}((cos(θ)(1 - \frac {1}{ρ^2}))\vec{g_θ} ) - (\frac{\partial}{\partial ρ}(sen (θ) (ρ + \frac {1}{ρ}))\vec{g_z})-(\frac{\partial}{\partial θ}((cos(θ)(1 - \frac {1}{ρ^2}))\vec{g_z} + (\frac{\partial}{\partial z}(sen (θ) (ρ + \frac {1}{ρ})\vec{g_ρ}) )=[/math] [math]\frac {1}{ρ}(-sen(θ) (1- \frac{1}{ρ^2}) \vec{g_z} + sen(θ) (1- \frac{1}{ρ^2}) \vec{g_z})= 0 [/math]

5 Divergencia de [math]\vec u[/math] nula

Que la divergencia de [math]\vec u[/math] sea nula se debe a que se trata de un fluido incompresible que implica que este tiene densidad constante en todo el conjunto y a lo largo del tiempo.

La divergencia mide la diferencia de flujo saliente y flujo entrante de un campo vectorial sobre la superficie que rodea un volumen control. Si la divergencia en un punto es positiva, se dice que el campo posee fuentes. Si la divergencia es negativa, se dice que tiene sumideros.

Vamos a comprobarlo analíticamente realizando los orrespondientes cálculo en coordenadas cilíndricas::

[math]\nabla·\vec{u}= \frac {1}{g^\frac{1}{2}}\frac{\partial }{\partial x^i}(g^\frac{1}{2} u^i)= [/math] [math]= \frac {1}{ρ}(\frac{\partial }{\partial ρ}(ρ u^ρ) + \frac{\partial }{\partial θ}(ρ u^θ) + \frac{\partial }{\partial z}(ρ u^z))=[/math] [math]\frac {1}{ρ}(\frac{\partial }{\partial ρ}((ρ - \frac{1}{ρ})cos(θ)) - \frac{\partial }{\partial θ}(1 + \frac{1}{ρ^2})sen(θ) + 0 )=[/math] [math]((1 + \frac{1}{ρ^2})cos(θ) - (1 + \frac{1}{ρ^2})cos(θ))= 0[/math]

6 Líneas de corriente

Las líneas de corriente del campo [math]\vec u[/math] son las líneas que son tangentes a [math]\vec u[/math] en cada punto. Para calcularlas primeramente hallamos el campo [math]\vec u[/math] que en cada punto es ortogonal a [math]\vec u[/math] , luego hacemos: [math]\vec{v}=\vec{k}\times\vec{u}[/math]

[math]\vec{v}=\vec{k}\times\vec{u}= {g^\frac{1}{2}}\begin{pmatrix} \vec{i} & \vec{j} & \vec{k} \\ 0 & 0 & 1 \\ u^1 & u^2 & u^3 \end{pmatrix}= ({u^1} \vec{j} - {u^2} \vec{i}) = (\frac{2xy}{(x^2+y^2 )^2})\vec{i} + (1 + \frac{y^2-x^2}{(x^2+y^2 )^2})\vec{j}[/math]

Este nuevo campo vectorial [math]\vec v[/math] es irrotacional, es decir que su rotaciones es 0; como demostramos a continuación ya que [math]\vec u[/math] es de divergencia nula.

[math]\nabla\times\vec{v}= \frac {1}{g^\frac{1}{2}}\begin{pmatrix} \vec{i} & \vec{j} & \vec{k} \\\frac{\partial }{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z} \\ v_1 & v_2 & v_3 \end{pmatrix}=(\frac{∂v_2}{∂x}-\frac{∂v_1}{∂y})\vec k = (\frac{∂u_1}{∂x}+\frac{∂u_2}{∂y}) \vec k =(\nabla·\vec u)\vec k = 0 \vec k =\vec 0[/math]

Y tiene un potencial escalar [math]ψ[/math] que se conoce como función de corriente de [math]\vec u[/math] y que vamos a calcular:

[math]\vec{v}=\vec{grad}ψ= \frac{∂ψ}{∂x}\vec{i} + \frac{∂ψ}{∂y}\vec{j} + \frac{∂ψ}{∂z}\vec{k} = (\frac{2xy}{(x^2+y^2 )^2})\vec{i} + (1 + \frac{y^2-x^2}{(x^2+y^2 )^2})\vec{j}[/math]

De esta fórmula anterior podemos sacar que:

[math]\frac{∂ψ}{∂x}= \frac{2xy}{(x^2+y^2 )^2}[/math]

[math]\frac{∂ψ}{∂y}= 1+ \frac{y^2-x^2}{(x^2+y^2 )^2}[/math]

Despejando [math]ψ[/math] de [math]\frac{∂ψ}{∂x}[/math] e Integrando obtenemos la función de [math]ψ(x,y)[/math]: [math]ψ=\int \frac{2xy}{(x^2+y^2 )^2} dx= -\frac{y}{x^2+y^2} + h(y)[/math]

A continuación derivamos la [math]ψ[/math] obtenida respecto "y" y la igualamos a la expresión que teníamos anterior de [math]\frac{∂ψ}{∂y}[/math] de manera que así podemos sacar [math]h(y)[/math] y completar la función [math]ψ[/math]

[math]\frac{∂ψ}{∂y}=\frac{y^2-x^2}{(x^2+y^2 )^2}+h'(y)= 1+ \frac{y^2-x^2}{(x^2+y^2 )^2}[/math]

Despejamos [math]h'(y)[/math] y sacamos integrándola [math]h(y)[/math]: [math]h'(y)=1[/math]

[math]h(y)=y[/math]

Por lo tanto el potencial escalar [math]ψ[/math] es igual a: [math]ψ=y-\frac{y}{x^2+y^2}[/math]

Dibujamos las líneas de corriente del campo [math]ψ_1[/math] y vemos que efectivamente son las líneas de corriente de [math]u_1[/math]:

u=linspace(1,5,100);
v=linspace(0,2*pi,100); 
[uu,vv]=meshgrid(u,v);
xx=uu.*cos(vv);       
yy=uu.*sin(vv);
fi=inline('y-y./(x.^2+y.^2)','x','y'); %representamos sobre el mallado
f=fi(xx,yy);
contour(xx,yy,f,30)
hold on
x=cos(v);
y=sin(v);
plot(x,y,'k','Linewidth',1)
axis([-4,4,-4,4])  % select region for drawing
axis equal
 
 
u=linspace(1,5,15);
v=linspace(0,2*pi,15); 
[uu,vv]=meshgrid(u,v); %mallado sobre el que representar
xx=uu.*cos(vv);       
yy=uu.*sin(vv);
fx=inline('1+(y.^2-x.^2)./((y.^2+x.^2).^2)','x','y');                 %componente x del vector 
fy=inline('-2*(x.*y)./((y.^2+x.^2).^2)','x','y');                 %compenente y del vector
u1=fx(xx,yy);
u2=fy(xx,yy);
quiver(xx,yy,u1,u2)     %dibujar el vector
view(2) 
hold off

Lineas de corriente del campo [math] \vec u [/math]

7 Estudio de la velocidad del fluido

A continuación se realiza el estudio de la velocidad del fluido en los puntos de la frontera S, que son aquellos en los que se particulariza la función del campo vectorial [math]\vec u[/math] para ρ=1.

Para determinar los valores máximo y mínimo de la velocidad consideramos su módulo:

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

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

Los valores máximo y mínimos se obtienen sustituyendo en la f(ρ,θ) para ρ=1:

[math]f(1,θ)=\sqrt((4(sen(θ))^2))=2|sen(θ)|[/math]

Por tanto, la velocidad será máxima cuando la función seno tome su máximo, es decir, cuando toma el valor 1:

Valores máximos y mínimos de la velocidad del fluido

[math]|\vec u|(\frac{ \pi}{ 2})= 2[/math]

[math]|\vec u| (\frac{ 3\pi}{ 2})= 2[/math]

En los puntos en coordenadas polares (1,π⁄2) y (1,3π⁄2), la velocidad máxima será igual a 2. Por su parte, la velocidad será mínima cuando la función seno tome su valor mínimo, es decir, cuando toma el valor 0: [math]|\vec u| (0)= 0 [/math]

[math]|\vec u| (\pi)= 0 [/math]

En los puntos en coordenadas polares (1,0) y (1,π), la velocidad mínima será igual a 0. Dichos puntos serán puntos de remanso, pues en ellos su velocidad es nula.

8 Estudio de la presión del fluído

Basándonos en la ecuación de Bernouilli:

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

y suponiendo que la densidad del fluido es constante [math]ρ = 2[/math]; vamos a calcular la presión del fluido dando un valor a la cte igual a 50: [math] p(ρ,⁡θ)=50-|\vec u|^2=50-((sin⁡θ (1+\frac{1}{ρ^2 })^2+(cos⁡θ(1-\frac{1}{ρ^2})^2 ) [/math] Para su representación expresamos la Presión en coordenadas cartesianas en lugar de cilíndricas: [math] p(x,⁡y)=50-|\vec u|^2=50-(1+ \frac{y^2-x^2}{(y^2+x^2)^2})^2+(\frac{-2xy}{(y^2+x^2)^2})^2) [/math] Y su representación en unos ejes de coordenadas planos respecto la Velocidad queda:

u=linspace(0,1,20);
v=linspace(0,1,20); 
fi=inline('(1+(y.^2-x.^2)./((y.^2+x.^2).^2)).^2+(-2*(x.*y)./((y.^2+x.^2).^2)).^2','x','y'); %Función Velocidad
f=fi(u,v);
p=50-f; %Función Presión
plot(f,p) %Representamos la función presión dependiendo de la velocidad
 
xlabel('velocidad')
ylabel('presion')
axis([0,50,0,50])

Gráfica de la ecuación de Bernouilli

Como podemos observar en el gráfico, a 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.

Si se consideran los puntos de mínima velocidad (v = 0), el primer término se anula y entonces la presión alcanza su valor máximo, que se corresponde con el valor de la constante. Por el contrario, si se toman los puntos de velocidad máxima se tendrá el menor valor de presión. Sabiendo que la velocidad máxima es 2 y la densidad del fluido es 2, el valor del primer término es 4 y se puede determinar que el valor de la constante tiene que ser mayor que 4 para que se obtenga un valor de presión positivo.

9 Trayectoria de una partícula del fluído

Supongamos que nuestra partícula va por una trayectoria horizontal paralela y próxima al eje del diámetro horizontal del objeto circular, por el hecho de estar en el fluido llevara una determinada velocidad y tendrá una determinada presión, al aproximarse al objeto empezara a curvar su trayectoria con el objetivo de rodear el objeto (en el caso de la partícula que choca justo con el centro del obstáculo llegará a detenerse); de manera que alcanzará su velocidad mínima al estar apunto de detenerse y irá acelerando, por lo que disminuirá su presión al sortear el obstáculo alcanzando su máxima velocidad al llegar a ser tengente a la circunferencia (corte con el diámetro vertical. Una vez alcanzado este punto empezará a disminuir su velocidad hasta alcanzar de nuevo el mínimo diametralmente opuesto al anteriormente descrito; de manera que su velocidad irá disminuyendo y su presión aumentando; hasta volver a una trayectoria horizontal paralela y próxima al eje con las condiciones de velocidad y presión previas a la acción del obstáculo.

10 Circulación a lo largo de la circunferencia unidad

La circulación del campo [math]\vec u[/math] a lo largo de la circunferencia unitaria centrada en el origen es nula como podemos observar a continuación debido a que esta circulación por el Teorema de Stokes se puede expresar como el flujo del rotacional a través de la superficie, y este rotacional como comprabamos anteriormente era 0:

[math] ∫ \vec {u} \vec {dr} =∫(\nabla\times\vec{u}) \vec {N} dS= ∫(\vec 0) \vec {N} dS= 0[/math]

Al ser el nula debido al teorema de Kutta-Joukowski, que establece que la fuerza que ejerce el fluido sobre el obstáculo (en nuestro caso el círculo unidad) es proporcional a la circulación y también a la densidad del fluido y la velocidad de este; resulta que le fluido no ejerce ninguna fuerza sobre el obstáculo. Esto es lo que se conoce como la paradoja de D'Alembert.

11 Repetición del análisis con [math]φ_1= (ρ+\frac{1}{ρ})cos(θ) + \frac{θ}{4π}[/math]

11.1 Velocidad de las partículas del fluido

La velocidad de las partículas del fluido viene dada por un campo vectorial [math]\vec u_1 [/math] cuya función potencial es [math]φ_1 (ρ,θ)=(ρ+\frac{1}{ρ})cos(θ) + \frac{θ}{4π}= φ_1 (x,y)=x+\frac{x}{x^2+y^2}+\frac{arctg(\frac{y}{x})}{4π}[/math] . Luego [math]\vec u_1[/math] es el gradiente de la función potencial: : [math]\vec u_1=∇φ_1[/math] Realizaremos los cálculos en cilíndricas y después lo pasaremos a cartesianas:

[math]∇φ_1 = \vec u=\frac{∂φ_1}{∂x^i} \vec g^i= [/math]

[math] =\frac{∂φ_1}{∂ρ} \vec g^ρ + \frac{∂φ_1}{∂θ} \vec g^θ + \frac{∂φ_1}{∂z} \vec g^z= [/math]

[math] =\frac{∂φ_1}{∂ρ} \vec g_ρ + \frac{1}{ρ^2} \frac{∂φ_1}{∂θ} \vec g_θ + \frac{∂φ_1}{∂z} \vec g_z = [/math]

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

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

Entonces las Componentes covariantes del vector quedan:

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

Y las Componentes contravariantes del vector :

[math] \vec u (ρ,θ)=cos(θ)(1 - \frac {1}{ρ^2}) \vec g_ρ - sen (θ) (\frac {1}{ρ} + \frac {1}{ρ^3})+ \frac{1}{4πρ^2} \vec g_θ [/math]

Pasadas las componentes a coordenadas cartesianas el vector [math]\vec u[/math] queda:

[math] \vec u (x,y)=(\frac{1+(y^2-x^2)}{(x^2+y^2 )^2} -\frac{y}{4πx^2(1+(\frac{y}{x})^2)} )\vec i +(\frac{1}{4πx(1+(\frac{y}{x})^2)}-\frac{2xy}{(x^2+y^2 )^2}) \vec j [/math]

Las representaciones de estos dos campos, el campo vectorial [math]\vec u_1 [/math] y el campo escalar φ_1, función potencial de [math]\vec u_1 [/math], son las siguientes:

%Representación del campo escalar φ_1
u=linspace(1,5,100);
v=linspace(0,2*pi,100); 
[uu,vv]=meshgrid(u,v);
xx=uu.*cos(vv);       
yy=uu.*sin(vv);
fi=inline('x+x./(x.^2+y.^2)+(1/(4*pi))*atan(y./x)','x','y'); %representamos sobre el mallado
f=fi(xx,yy);
contour(xx,yy,f,30)
axis([-4,4,-4,4])      % select region for drawing
hold on
x=cos(v);
y=sin(v);
plot(x,y,'k','Linewidth',1)
axis equal
 
%Representación del campo vectorial u_1
u=linspace(1,5,15);
v=linspace(0,2*pi,15); 
[uu,vv]=meshgrid(u,v); %mallado sobre el que representar
xx=uu.*cos(vv);       
yy=uu.*sin(vv);
fx=inline('1+(y.^2-x.^2)./((y.^2+x.^2).^2)-y./(4*pi*(x.^2).*(1+(y./x).^2))','x','y');                 %componente x del vector 
fy=inline('1./(4*pi*x.*(1+(y./x).^2))-2*(x.*y)./((y.^2+x.^2).^2)','x','y');                 %compenente y del vector
u1=fx(xx,yy);
u2=fy(xx,yy);
quiver(xx,yy,u1,u2)     %dibujar el vector
 
hold off

Representación de la función potencial y del campo escalar [math] \vec u_1[/math]

Como podemos observar, el campo vectorial [math]\vec u_1 [/math] de las velocidades de las partículas es ortogonal a las curvas de nivel de [math]φ_1[/math] ; porque al ser el gradiente de la función potencial representa la línea de máxima pendiente del campo escalar [math]φ_1[/math].

11.2 Rotacional de [math]\vec u_1[/math] nulo

El rotacional de un campo nos indica la tendencia de este a inducir rotación alrededor de un punto.

El rotacional de nuestro campo [math]\vec u_1[/math] es nulo ya que todo campo potencial, es decir, expresable como el gradiente de un potencial escalar es irrotacional.

Vamos a comprobarlo entonces analíticamente realizando los correspondientes cálculo en coordenadas cilíndricas:

[math]\nabla\times\vec{u_1}= \frac {1}{g^\frac{1}{2}}\begin{pmatrix} \vec{g_ρ} & \vec{g_θ} & \vec{g_z} \\ \frac{ \partial}{\partial ρ} & \frac{\partial }{\partial θ} & \frac{\partial }{\partial z}\\ u1_1 & u1_2 & u1_3 \end{pmatrix}= \frac {1}{ρ} \begin{pmatrix} \vec{g_ρ} & \vec{g_θ} & \vec{g_z} \\ \frac{ \partial}{\partial ρ} & \frac{\partial }{\partial θ} & \frac{\partial }{\partial z}\\ cos(θ)(1 - \frac {1}{ρ^2})& - (sen (θ) (ρ + \frac {1}{ρ})) + \frac{1}{4π} & 0 \end{pmatrix}=[/math]

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

11.3 Divergencia de [math]\vec u_1[/math] nula

Que la divergencia de [math]\vec u_1[/math] sea nula se debe a que se trata de un fluido incompresible que implica que este tiene densidad constante en todo el conjunto y a lo largo del tiempo.

La divergencia mide la diferencia de flujo saliente y flujo entrante de un campo vectorial sobre la superficie que rodea un volumen control. Si la divergencia en un punto es positiva, se dice que el campo posee fuentes. Si la divergencia es negativa, se dice que tiene sumideros.

Vamos a comprobarlo analíticamente realizando los orrespondientes cálculo en coordenadas cilíndricas:

[math]\nabla·\vec{u_1}= \frac {1}{g^\frac{1}{2}}\frac{\partial }{\partial x^i}(g^\frac{1}{2} u1^i)= [/math]

[math]= \frac {1}{ρ}(\frac{\partial }{\partial ρ}(ρ u1^ρ) + \frac{\partial }{\partial θ}(ρ u1^θ) + \frac{\partial }{\partial z}(ρ u1^z))=[/math]

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

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

11.4 Líneas de corriente

Las líneas de corriente del campo [math]\vec u_1[/math] son las líneas que son tangentes a [math]\vec u_1[/math] en cada punto. Para calcularlas primeramente hallamos el campo [math]\vec v_1[/math] que en cada punto es ortogonal a [math]\vec u_1[/math] , luego hacemos: [math]\vec{v_1}=\vec{k}\times\vec{u_1}[/math]

[math]\vec{v_1}=\vec{k}\times\vec{u_1}= {g^\frac{1}{2}}\begin{pmatrix} \vec{i} & \vec{j} & \vec{k} \\ 0 & 0 & 1 \\ u1^1 & u1^2 & u1^3 \end{pmatrix}= 1 (\vec{u1^1} \vec{j} - \vec{u1^2} \vec{i}) = (\frac{2xy}{(x^2+y^2 )^2}-\frac{1}{4πx(1+(\frac{y}{x})^2)})\vec{i} + (1 + \frac{y^2-x^2}{(x^2+y^2 )^2}-\frac{y}{4πx^2(1+(\frac{y}{x})^2)})\vec{j}[/math]

Este nuevo campo vectorial [math]\vec v_1[/math] es irrotacional, es decir que su rotaciones es 0; como demostramos a continuación ya que [math]\vec u_1[/math] es de divergencia nula.

[math]\nabla\times\vec{v_1}= \frac {1}{g^\frac{1}{2}}\begin{pmatrix} \vec{i} & \vec{j} & \vec{k} \\\frac{\partial }{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z} \\ v1_1 & v1_2 & v1_3 \end{pmatrix}=(\frac{∂v1_2}{∂x}-\frac{∂v1_1}{∂y})\vec k = (\frac{∂u1_1}{∂x}+\frac{∂u1_2}{∂y}) \vec k =(\nabla·\vec u)\vec k = 0 \vec k =\vec 0[/math]

Y tiene un potencial escalar [math]ψ_1[/math] que se conoce como función de corriente de [math]\vec u_1[/math] y que vamos a calcular:

[math]\vec{v_1}=\vec{grad}ψ_1= \frac{∂ψ_1}{∂x}\vec{i} + \frac{∂ψ_1}{∂y}\vec{j} + \frac{∂ψ_1}{∂z}\vec{k} = (\frac{2xy}{(x^2+y^2 )^2}-\frac{1}{4πx(1+(\frac{y}{x})^2)})\vec{i} + (1 + \frac{y^2-x^2}{(x^2+y^2 )^2}-\frac{y}{4πx^2(1+(\frac{y}{x})^2)})\vec{j}[/math]

De esta fórmula anterior podemos sacar que:

[math]\frac{∂ψ_1}{∂x}= (\frac{2xy}{(x^2+y^2 )^2}-\frac{1}{4πx(1+(\frac{y}{x})^2)})[/math]

[math]\frac{∂ψ_1}{∂y}= (1 + \frac{y^2-x^2}{(x^2+y^2 )^2}-\frac{y}{4πx^2(1+(\frac{y}{x})^2)})[/math]

Despejando [math]ψ_1[/math] de [math]\frac{∂ψ_1}{∂x}[/math] e Integrando obtenemos la función de [math]ψ_1(x,y)[/math]: [math]ψ_1=\int \frac{2xy}{(x^2+y^2 )^2}-\frac{1}{4πx(1+(\frac{y}{x})^2)} dx= -\frac{y}{x^2+y^2}-\frac{Ln(x^2+y^2)}{8π} + h(y)[/math]

A continuación derivamos la [math]ψ_1[/math] obtenida respecto "y" y la igualamos a la expresión que teníamos anterior de [math]\frac{∂ψ_1}{∂y}[/math] de manera que así podemos sacar [math]h(y)[/math] y completar la función [math]ψ_1[/math]

[math]\frac{∂ψ_1}{∂y}=\frac{y^2-x^2}{(x^2+y^2 )^2}-\frac{2y}{8π(x^2+y^2)}+h'(y)= (1 + \frac{y^2-x^2}{(x^2+y^2 )^2}-\frac{y}{4πx^2(1+(\frac{y}{x})^2)})[/math]

Despejamos [math]h'(y)[/math] y sacamos integrándola [math]h(y)[/math]: [math]h'(y)=1[/math]

[math]h(y)=y[/math]

Por lo tanto el potencial escalar [math]ψ_1[/math] es igual a: [math]ψ_1=y-\frac{y}{x^2+y^2}-\frac{Ln(x^2+y^2)}{8π}[/math]

Dibujamos las líneas de nivel del campo [math]ψ_1[/math] y vemos que efectivamente son las líneas de corriente de [math]u_1[/math]:

u=linspace(1,5,100);
v=linspace(0,2*pi,100); 
[uu,vv]=meshgrid(u,v);
xx=uu.*cos(vv);       
yy=uu.*sin(vv);
fi=inline('y-y./(x.^2+y.^2)-(1/(8*pi))*log(x.^2+y.^2)','x','y'); %representamos sobre el mallado
f=fi(xx,yy);
contour(xx,yy,f,30)
hold on
x=cos(v);
y=sin(v);
plot(x,y,'k','Linewidth',1)
axis([-4,4,-4,4])  % select region for drawing
axis equal 

 
u=linspace(1,5,15);
v=linspace(0,2*pi,15); 
[uu,vv]=meshgrid(u,v); %mallado sobre el que representar
xx=uu.*cos(vv);       
yy=uu.*sin(vv);
fx=inline('1+(y.^2-x.^2)./((y.^2+x.^2).^2)-y./(4*pi*(x.^2).*(1+(y./x).^2))','x','y');                 %componente x del vector 
fy=inline('1./(4*pi*x.*(1+(y./x).^2))-2*(x.*y)./((y.^2+x.^2).^2)','x','y');                 %compenente y del vector
u1=fx(xx,yy);
u2=fy(xx,yy);
quiver(xx,yy,u1,u2)     %dibujamos el vector
 
hold off

Líneas de corriente de [math]\vec u_1[/math]