Visualización de campos escalares y vectoriales en fluidos. Grupo 8C

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en fluidos. Grupo 8C
Asignatura	Teoría de Campos
Curso	2021-22
Autores	Celia Vives Mora Jesús Gil Gutiérrez Pablo Martín Mulero
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

En el siguiente trabajo trataremos la visualización de una serie de campos escalares y vectoriales sobre la superficie de un canal rectangular. Para la representación gráfica de estos campos recurriremos a Octave.

1 Mallado del canal

En primer lugar, vamos a visualizar nuestro canal y comprender en dónde vamos a trabajar. Para ello, definiremos el mallado del canal sobre el que trabajaremos. Este canal posee una longitud y un ancho de 8 y 2 metros respectivamente. El código Octave es el que se adjunta a continuación:

x=0:0.1:8; %Vector X.
y=-1:0.1:1; %Vector Y.
[xx,yy]=meshgrid(x,y); %Mallado XY.
hold on
mesh(xx,yy,0.*xx); %Representamos el canal.
axis([0,8,-2,2]); %Rango de los ejes
axis equal %Los ejes X e Y iguales
view(2);
title ('Mallado del canal');
hold off

Gráficamente, obtenemos lo siguiente:

Mallado del canal

2 Ecuación de Navier-Stokes

Se nos dice que la velocidad del fluido en el canal viene dado por el siguiente campo vectorial: [math]\ \vec u (x,y)=(y+1)(1-y)((p_1-p_2)/2 {\mu})\vec i [/math] ; mientras que la presión viene dado por el siguiente campo escalar: [math]\ p(x,y)=p_{\ 1}+(p_{\ 2}-p_{\ 1})(x-1)[/math]. Sean [math]\ p_{\ 1} [/math] la presión en los puntos [math]\ x=1 [/math], [math]\ p_{\ 2} [/math] la presión en los puntos [math]\ x=2 [/math] y[math]\ {\mu} [/math] el coeficiente de viscosidad del fluido.

Las ecuaciones de Navier-Stokes describen el comportamiento de los fluidos "newtonianos", entre otros la conservación del movimiento y la conservación de la masa. Aplicando los principios de conservación de la mecánica y la termodinámica se obtienen estas fórmulas. En nuestro caso, la ecuación [math]\ (\vec u · ∇) \vec u +∇p={\mu}∆\vec u [/math] representa el principio de conservación del momento lineal. Vamos a comprobarlo:

[math]\ (\vec u·∇) \vec u= u_j \frac {\partial u_i}{\partial x_j} \vec e_i= u_1 \frac {\partial u_i}{\partial x_1} \vec e_i+ u_2 \frac {\partial u_i}{\partial x_2} \vec e_i=u_1 \frac {\partial u_1}{\partial x_1} \vec e_1+ u_1 \frac {\partial u_2}{\partial x_1} \vec e_2+ u_2 \frac {\partial u_1}{\partial x_2} \vec e_1+ u_2 \frac {\partial u_2}{\partial x_2} \vec e_2 [/math]; sean en nuestro caso [math]\ u_1=(y+1)(1-y)( \frac {p_1-p_2}{2 {\mu}}); u_2=0 [/math].

Sustituimos y tenemos lo siguiente: [math]\ u_1 \frac {\partial u_1}{\partial x_1} \vec e_1+ u_1 \frac {\partial 0}{\partial x_1} \vec e_2+ 0 \frac {\partial u_1}{\partial x_2} \vec e_1+ 0 \frac {\partial 0}{\partial x_2} \vec e_2= (y+1)(1-y)( \frac {p_1-p_2}{2 {\mu}})( \frac {\partial ((y+1)(1-y)( \frac {p_1-p_2}{2 {\mu}}))}{\partial x_1}) \vec e_1= (y+1)(1-y)( \frac {p_1-p_2}{2 {\mu}})(0) \vec e_1=0 [/math].

Vayamos ahora con el segundo término de la ecuación, que es el gradiente de p: [math]\ ∇p= \frac{\partial p}{\partial x_i} \vec e_i= \frac{\partial p}{\partial x} \vec i + \frac {\partial p}{\partial y} \vec j= \frac{\partial (p_1+(p_2-p_1)(x-1))}{\partial x} \vec i +0 \vec j= (p_2-p_1) \vec i [/math]

Nos queda el último término de la ecuación: [math]\ ∆ \vec u=∆u_i \vec e_i=∆u_1 \vec e_1+∆u_2 \vec e_2=( \frac {\partial^2 u_1}{\partial x_1^2}+ \frac {\partial^2 u_1}{\partial x_2^2}) \vec e_1+( \frac {\partial^2 u_2}{\partial x_1^2}+ \frac {\partial^2 u_2}{\partial x_2^2}) \vec e_2= (0+ \frac{\partial^2 ((y+1)(1-y)((p_1-p_2)/2 {\mu}))}{\partial x_2^2}) \vec e_1+ (0+0) \vec e_2= - \frac {p_1-p_2}{{\mu}} \vec e_1=- \frac {p_1-p_2}{{\mu}} \vec i [/math]; Llegamos entonces a [math]\ ∆ \vec u= \frac {p_2-p_1}{{\mu}} \vec i [/math]

Volvamos a la ecuación de Navier-Stokes: [math]\ (\vec u · ∇) \vec u +∇p={\mu}∆\vec u [/math], si sustituimos todo lo calculado anteriormente, llegamos a lo siguiente: [math]\ 0+(p_2-p_1) \vec i= {\mu}- \frac {p_1-p_2}{{\mu}} \vec i= -p_1+p_2 \vec i[/math]. Se cumple la Ecuación de Navier-Stokes Estacionaria. Es un fluido Newtoniano.

Vamos a verificar la condición de incompresibilidad [math]\ ∇·\vec u =0 [/math], veámoslo:

Tenemos [math]\ \vec u (x,y)=\vec u (x_{\ 1},x_{\ 2})=(x_{\ 2}+1)(1-x_{\ 2})(p_{\ 1}-p_{\ 2}/2 {\mu})\vec i [/math], sea la divergencia por definición [math]\ ∇·\vec u = \frac{\partial u_{\ 1}}{\partial x_{\ 1}}+\frac{\partial u_{\ 2}}{\partial x_{\ 2}}[/math]; operando: [math]\ ∇·\vec u = \frac{\partial(x_{\ 2}+1)(1-x_{\ 2})(p_{\ 1}-p_{\ 2}/2 {\mu})}{\partial x_{\ 1}}+\frac{\partial 0}{\partial x_{\ 2}}=0+0=0[/math]. Luego se cumple la condición. ¿Que quiere decir que la divergencia sea nula? Pues que el fluido ni se contrae ni se expande, osea, es incompresible.

La velocidad del fluido en las paredes es la velocidad en [math]\ y=-1 [/math] e [math]\ y=1 [/math]:

Velocidad en [math]\ y=-1 [/math] es: [math]\ \vec u (x_{\ 1},-1)=(0)(2)(\frac{p_{\ 1}-p_{\ 2}}{2 {\mu}})\vec i =0 [/math]

Velocidad en [math]\ y=1 [/math] es: [math]\ \vec u (x_{\ 1},1)=(2)(0)(\frac{p_{\ 1}-p_{\ 2}}{2 {\mu}})\vec i =0 [/math]

Osea, la velocidad del fluido en las paredes del canal es 0, esto se debe al rozamiento con estas.

3 Campo de presiones y velocidades

Recordemos que la presión venía dada por[math]\ p(x,y)=p_{\ 1}+(p_{\ 2}-p_{\ 1})(x-1)[/math] y la velocidad por [math]\ \vec u (x,y)=(y+1)(1-y)(\frac{p_1-p_2}{2{\mu}})\vec i [/math], que cumplían la ecuación de Navier-Stokes estacionaria. Vamos a dar ahora unos valores a [math]\ p_{\ 1} [/math],[math]\ p_{\ 2} [/math] y[math]\ {\mu} [/math]. Como son [math]\ p_{\ 1}=2 [/math],[math]\ p_{\ 2}=1 [/math] y [math]\ {\mu}=1 [/math]. Dibujemos el campo de presiones y velocidades asociados a estos valores:

Campo de presiones

x=0:0.1:8;
y=-1:0.1:1
[xx,yy]=meshgrid(x,y);
p=2+(1-2)*(xx-1);
hold on
surf(xx,yy,p);
axis([0,8,-2,2]);
axis equal
colorbar
title('Campo de presiones')
hold off

Podemos apreciar que cuanto más nos alejamos del origen las presiones irán descendiendo.

Campo de velocidades

x=0:0.4:8; %Separamos la longitud del paso para que se vea mejor el campo vectorial
y=-1:0.1:1;
[xx,yy]=meshgrid(x,y);
vx=(yy+1).*(1-yy).*(2-1)./(2*1);
vy=0.*yy;
hold on
quiver(xx,yy,vx,vy);
axis([0,8,-2,2]);
axis equal
title('Campo de velocidades')
hold off

Con esta gráfica, vemos lo que ya habíamos calculado en el apartado anterior: La velocidad del fluido en las paredes es nula.

4 Puntos con mayor módulo de velocidad

Calculemos los puntos con mayor módulo de velocidad:

[math]\ |\vec u|= \frac {1}{2}-y^2 \frac {1}{2} [/math]

Si [math]\ y=0; \vec u=1 [/math]

Si [math]\ y=-1; \vec u=0 [/math]

Si [math]\ y=1; \vec u=0 [/math]

Si [math]\ y=0'5; \vec u=0'75 [/math]

Luego la velocidad es máxima en el centro del canal, es decir, en [math]\ y=0 [/math].

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

Hemos dibujado el campo de velocidades en el apartado anterior, pero ¿Cómo son las líneas de corriente?. Una línea de corriente es el camino que seguirá una partícula a lo largo del fluido, se trata del vector tangente al campo de velocidades. Veamos que [math] \vec v [/math] es irrotacional: Para ello calculamos el campo [math]\ \vec v [/math] que en cada punto es ortogonal a [math]\ \vec u [/math]. Tomaremos [math]\ \vec v= \vec k× \vec u [/math].

[math]\ \vec v= \vec k× \vec u= \vec k×((y+1)(1-y)( \frac {p_1-p_2}{2 {\mu}})) \vec i=(-y^2+1( \frac {p_1-p_2}{2 {\mu}})) \vec j [/math]

Sustituyendo [math]\ p_1=2, p_2=1, {\mu}=1 [/math], tenemos entonces [math]\ \vec v=(-y^2+1) \vec j [/math]

¿Es irrotacional? Veámoslo: [math]\ ∇× \vec v=\left|\begin{matrix} \vec i & \vec j & \vec k \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ 0 & -y^2+1 & 0 \end{matrix}\right|=0 \vec i+0 \vec j+\frac {\partial (-y^2+1)}{\partial x} \vec k-(\frac {\partial (-y^2+1)}{\partial z} \vec i+0 \vec j+0 \vec k)=0 [/math]. Luego es irrotacional. Es irrotacional ya que [math]\ \vec u [/math] tiene divergencia nula.

Vamos a calcular [math]\ ψ: \vec v=∇ψ=\frac{\partial ψ}{\partial x}\vec e_i+\frac{\partial ψ}{\partial y}\vec e_j=0+\frac{\partial ψ}{\partial y} [/math]; Vamos a integrar: [math]\ \int \frac{p_1-p_2}{2{\mu}}dy-\int \frac{y^2(p_1-p_2)}{2 {\mu}}dy [/math]; [math]\ ψ=(\frac{y(p_1-p_2)}{2{\mu}}-\frac{y^3}{1}·\frac{(p_1-p_2)}{6{\mu}}); ψ=\frac {3y(p_1-p_2)-y^3(p_1-p_2)}{6{\mu}}[/math]

Líneas de corriente

x=0:0.25:8;
y=-1:0.2:1;
[X,Y]=meshgrid(x,y);
figure(4)
hold on
psi=Y/2-Y.^3/2;
contour(X,Y,psi,10,'r')
grid
axis([0,8,-2,2]);
axis equal;
title('Líneas de corriente')
xlabel('Eje x')
ylabel('Eje y')
view(2)
colorbar
quiver3(X,Y,0*X,1/2*X./X-Y.^2/2,0*X,0*X,'k');
hold off

En la gráfica de la derecha podemos ver el campo de velocidades [math]\ \vec u [/math] así como las líneas de corriente en rojo. Nótese que las líneas de corriente son rectas que siguen la misma dirección que el campo de velocidades. Esto se debe a que el campo se mueve en una única dirección.

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

[math]\ ∇× \vec u\left|\begin{matrix} \vec i & \vec j & \vec k \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ u_1 & u_2 & u_3 \end{matrix}\right| = \left|\begin{matrix} \vec i & \vec j & \vec k \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ \frac{(y+1)(y-1)(p_1-p_2)}{2{\mu}} & 0 & 0 \end{matrix}\right| =[/math] [math]\ 0 \vec j-\frac{\partial}{\partial y}(\frac{(y+1)(y-1)(p_1-p_2)}{2{\mu}})\vec k=\frac{(-2y)(p_1-p_2)}{2{\mu}}\vec k;|∇× \vec u|=\frac{-(p_1-p_2)}{{\mu}}=0[/math]

Como no existen ni máximos ni mínimos absolutos, tendremos que ver si las paredes del canal son máximos o mínimos relativos. Estamos trabajando en el intervalo y=[-1,1]. [math]\ y=1; |∇× \vec u|=|\frac{(-2·(-1))(p_1-p_2)}{2{\mu}}|=|\frac{p_1-p_2}{{\mu}}|;[/math] si [math]\ p_1=2 , p_2=1 , {\mu}=1[/math], entonces llegamos a lo siguiente: [math]\ |∇× \vec u|=|\frac{p_1-p_2}{{\mu}}|=|\frac{2-1}{1}|=|1|=1[/math]

[math]\ y=-1; |∇× \vec u|=|\frac{(-2·(1))(p_1-p_2)}{2{\mu}}|=|\frac{-(p_1-p_2)}{{\mu}}|;[/math] si [math]\ p_1=2 , p_2=1 , {\mu}=1[/math], entonces llegamos a lo siguiente: [math]\ |∇× \vec u|=|\frac{-(p_1-p_2)}{{\mu}}|=|\frac{-(2-1)}{1}|=|-1|=1[/math]

Rotacional de la velocidad

Con ello, concluimos que nuestros puntos dentro del canal con mayor rotacional son las paredes. Veámoslo de forma gráfica a continuación en Octave:

x=0:0.1:8;
y=-1:0.1:1;
[xx,yy]=meshgrid(x,y);
rot=abs(yy);
hold on
surf(xx,yy,rot);
colorbar;
view(2);
axis([0,8,-2,2]);
axis equal
title('Rotacional de la velocidad')
hold off

Se ve en la gráfica de la derecha que las paredes son los puntos de mayor rotacional, donde el valor es 1. Esto tiene una explicación: El rotacional mide lo que tiende a rotar un objeto/partícula debido al efecto del campo de velocidades; en el centro del canal -donde la velocidad es máxima- el descenso a izquierda o derecha (nos es indiferente, el campo es simétrico respecto al eje [math]\ y=0 [/math]) es a nivel infinitesimal cercano a 0. Sin embargo, el campo de velocidades sigue una progresión parabólica, que hace que el descenso sea cada vez más acusado según nos acercamos a los laterales, donde es máximo. Como es en los laterales ese descenso infinitesimal máximo, una partícula tenderá a rotar más en esa zona.

Representación esquemática del canal y el campo de velocidades

7 Caudal circulante por el centro del canal

El porcentaje de caudal circulante por el centro del canal viene dado por la siguiente relación: [math]\ {\%} Q_{central}= \frac{Q_{centro}}{Q_{total}}·100 [/math]. Muy bien, necesitamos saber entonces [math]\ Q_{centro} [/math] y [math]\ Q_{total} [/math]. Calculémoslo, para hacerlo vamos a plantearlo mediante integrales de superficie: Tomamos el caudal [math]\ \vec u(x,y)=(y+1)(1-y)(\frac{p_1-p_2}{2 {\mu}}) [/math], la superficie a atravesar será la sección del canal. Osea: [math]\ (y,z)∈(-1,1)×(0, \frac {1}{2}) [/math]; si parametrizamos esta sección tenemos lo siguiente: [math]\ y=u [/math], [math]\ z=v [/math], con los parámetros moviéndose en [math]\ (u,v)∈(-1,1)×(0, \frac {1}{2}) [/math]. Hecha esta parametrización, obtenemos el vector que la define: [math]\ \vec s(u,v)=u \vec j + v \vec k [/math]. Para simplificar un poco los cálculos nombraremos [math]\ q [/math] al término [math]\ \frac {p_1-p_2}{2 {\mu}} [/math] del campo vectorial.

La integral de superficie será la siguiente: [math]\ Q_{total}= \int \int \vec u( \vec s(u,v))·( \frac{\partial \vec s}{\partial u}× \frac {\partial \vec s}{\partial v})dvdu=\int_{u=-1}^{u=1} \int_{v=0}^{v=1/2}q(-u^2+1) \vec i ·(\vec j × \vec k) = q \int_{u=-1}^{u=1} \int_{v=0}^{v=1/2}q(-u^2+1) \vec i · \vec i= [/math] [math]\ q \int_{-1}^{1} (-u^2+1) \frac {1}{2} du= \frac {q}{2}(-[ \frac {u^3}{3}]_{u=-1}^{u=1}+[u]_{u=-1}^{u=1})[/math]. Terminando de operar, llegamos a lo siguiente: [math]\ \frac {q}{2}( \frac {-2}{3}+2)= \frac {2}{3}q[/math] Hacemos lo mismo pero con la parte central del canal, luego lo único que cambia son los límites de integración, en este caso [math]\ u∈(-0.5,0.5)[/math]. Con ello, el flujo del campo vectorial en la parte central es [math]\ Q_{centro}= \frac {11}{24}q [/math].

Calculemos el porcentaje de caudal por el centro del canal [math]\ {\%} Q_{central}= \frac{Q_{centro}}{Q_{total}}·100= \frac{11/24·q}{2/3·q}·100=68'75 {\%}. [/math].

Nótese que por la mitad central del canal no circula la mitad del caudal total, sino que es un porcentaje realmente elevado. Esto es debido a que la velocidad del fluido en el canal no es uniforme, si no que es máxima en el centro, descendiendo a partir de ahí siguiendo una progresión parabólica hasta las paredes, donde la velocidad es nula.

8 Campo de temperaturas

Se nos dice que la temperatura viene dada por el siguiente campo escalar en coordenadas polares: [math]\ T({\rho},{\theta})=1+{\rho}^2 sen^2 {\theta}e^{-({\rho-1/2})^2} [/math]. Lo primero de todo, pasemos el campo a coordenadas cartesianas, aplicando los siguientes cambios de variable: [math]\ {\rho}={\sqrt {x^2+y^2}} [/math] y [math] {\theta}=arctan(\frac {y}{x+0.0000001}) [/math]. El sumar 0.0000001 en el denominador es para evitar que nos de la indeterminación partido de 0 el programa. Teniendo en cuenta este cambio de variable, vamos a dibujar el campo de temperaturas y sus curvas de nivel, siendo el código Octave es el siguiente:

Temperatura del fluido

x=0:0.1:8;
y=-1:0.1:1;
[X,Y]=meshgrid(x,y);
rho=sqrt(X.^2+Y.^2);
hold on
T=X./X+Y.^2.*[exp(-(rho-1/2*X./X).^2)];
surf(X,Y,T);
view
axis([0,8,-2,2]);
axis equal
colorbar
title('Campo de temperaturas')
hold off

Para obtener las curvas de nivel de la temperatura, haremos lo siguiente:

Curvas de nivel de la temperatura

x=0:1/30:8;
y=-1:1/20:1;
[X,Y]=meshgrid(x,y);
rho=sqrt(X.^2+Y.^2);
hold on
contour(X,Y,T,10);
axis equal;
axis([0,8,-2,2]);
view(2);
title('Curvas de nivel de temperatura')
colorbar
hold off

Gráficamente, podemos ver que las zonas de mayor temperatura están en torno a [math]\ x=0[/math], junto a las paredes del canal.

Por último vamos a calcular el gradiente de la temperatura, de esa forma veremos en que dirección tiene un mayor crecimiento. Sea la temperatura dada en coordenadas polares [math]\ T({\rho},{\theta})=1+{\rho}^2 sen^2 {\theta}e^{-({\rho-1/2})^2} [/math]. Luego el gradiente en coordenadas polares, es por definición lo siguiente: [math]\ ∇T=\frac {\partial T({\rho},{\theta})}{\partial {\rho}} \vec i+ \frac {\partial T({\rho},{\theta})}{\partial {\theta}} \vec j [/math]; Aplicando la fórmula: [math]\ ∇T= \frac {\partial (1+{\rho}^2 sen^2 {\theta}e^{-({\rho-1/2})^2})}{\partial {\rho}} \vec i + \frac{\partial (1+{\rho}^2 sen^2 {\theta}e^{-({\rho-1/2})^2})}{\partial {\theta}} \vec j=(2 {\rho}sen^2{\theta}e^{-({\rho}^2- \frac {1}{2})^2}) \vec i-(4{\rho({\rho}^2- \frac {1}{2})e^{-({\rho}^2- \frac {1}{2})^2}}) \vec j[/math]

Si dibujamos las curvas de nivel y el gradiente, obtenemos lo siguiente:

Gradiente de la temperatura

x=0:0.1:8;
y=-1:0.1:1;
[X,Y]=meshgrid(x,y);
rho=sqrt(X.^2+Y.^2);
figure(2)
hold on
quiver(X,Y,-2*X.*Y.^2.*(rho-1/2*X./X).*exp(-(rho-1/2*X./X).^2)./rho,2*Y.*exp(-(rho-1/2*X./X).^2)-2*Y.^3.*(rho-1/2*X./X).*exp(-(rho-1/2*X./X).^2)./rho);
contour(X,Y,T,10)
axis equal
grid
view(2)
axis equal
axis([0,8,-2,2]);
colorbar
title('Gradiente temperatura y curvas de nivel')
hold off

Nótese que el gradiente es ortogonal a las curvas de nivel siempre, veámoslo con mayor ampliación:

Gradiente de la temperatura y curvas de nivel

El gradiente indica la dirección de máximo crecimiento de una función, mientras que las curvas de nivel son los puntos donde la función escalar tiene un valor constante, por ello siempre van a ser ortogonales.