Visualización de campos escalares y vectoriales en fluidos (Grupo C20)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Visualización de campos escalares y vectoriales en fluidos. (Grupo C20) |
| Asignatura | Teoría de Campos |
| Curso | 2014-15 |
| Autores | Javier Ruiz de Galarreta López, Luis Moreno Fernandez de Soria, Eduardo Moreda Meleiro, Jose Antonio Martínez Montalvo, Alberto Rodríguez Ruiz, Natalia Opie Dávila |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
1 Introducción
1.1 Objetivos y metodología
Este trabajo tiene como objetivo el estudio de los campos físicos escalares y vectoriales que definen un fluido que discurre por un canal. Este estudio se ha realizado en dos dimensiones, y para la visualización de los campos y el cálculo numérico de algunas integrales se han empleado los programas Matlab y Octave UPM.
El campo que vamos a tratar viene definido por los siguientes campos:
Campo vectorial de velocidades: [math]\vec{u}(x,y)=u_1(x,y)\vec{i}+u_2(x,y)\vec{j}=y\cdot(1-y)\frac{p_1-p_2}{2μ}\vec{i}[/math]
Campo escalar de presiones: [math]p(x,y)=p_1+(p_2-p_1)(x-1)[/math]
Campo escalar de temperatura: [math]T(x,y)=e^{-(x-1)^2+y^2}[/math]
Donde [math]μ[/math] es el coeficiente de viscosidad del fluido, [math]p_1[/math] es la presión en los puntos x = 1 y [math]p_2[/math] es la presión en los puntos x = 2
1.2 Representación del espacio ocupado por el fluido
El fluido objeto del estudio se halla ocupando el espacio delimitado por el rectángulo [0,4] x [-1,2].
Para representar gráficamente dicho dominio, crearemos el siguiente programa en MATLAB:
% Representación gráfica del dominio ocupado por un fluido
% mediante un mallado que abarca los puntos interiores del
% rectangulo [0; 4]x[0; 1], con los ejes fijados en la
% region [0; 4][-1; 2].
hold off
% Mallado:
[X,Y] = meshgrid([0:0.1:4],[0:0.1:1]);
mesh(X,Y,0*X)
% Ejes, región y rectas superior e inferior:
grid on
axis([0,4,-1,2])
axis equal
xlabel('X')
ylabel('Y')
hold on
x = [0,4];
y = [0,0];
plot(x,y,'-k','linewidth',2)
hold on
x = [0,4];
y = [1,1];
plot(x,y,'-k','linewidth',2)
Las ecuaciones de Navier-Stokes son un conjunto de ecuaciones en derivadas parciales no lineales que describen el movimiento de un fluido. Estas ecuaciones gobiernan la atmósfera terrestre, las corrientes oceánicas y el flujo alrededor de vehículos o proyectiles y, en general, cualquier fenómeno en el que se involucren fluidos newtonianos. Reciben su nombre de Claude-Louis Navier (1785-1836) y George Gabriel Stokes (1819-1903).
En este caso vamos a trabajar con la ecuación de Navier Stokes estacionaria, cuya expresión es la que sigue:
Donde [math]\vec{u}(x,y)[/math] y [math]p(x,y)[/math] son los campos de velocidades y presiones del fluido descritos en la introducción, y [math]μ[/math] es el coeficiente de viscosidad del fluido.
La misma ecuación empleando índices será:
Operando se tiene:
- El gradiente de [math]\vec{u}(x,y)[/math] :
[math]((\nabla\vec{u})_{ij})= \begin{pmatrix} \frac{\partial u_1}{\partial x} & \cfrac{\partial u_1}{\partial y} \\ \frac{\partial u_2}{\partial x} & \cfrac{\partial u_2}{\partial y} \\ \end{pmatrix} = \begin{pmatrix} 0 & (1-2y)\frac{p_1-p_2}{2μ} \\ 0 & 0 \\ \end{pmatrix} [/math]
- El gradiente de [math]p(x,y)[/math] :
[math]\begin{pmatrix}
(\nabla\ p)_1 \\
(\nabla\ p)_2 \end{pmatrix} = \begin{pmatrix}
\frac{\partial p}{\partial x} \\
\frac{\partial p}{\partial y} \end{pmatrix} = \begin{pmatrix}
p_2-p_1 \\
0 \end{pmatrix} [/math]
- El laplaciano de [math]\vec{u}(x,y)[/math]:
[math]\Delta\vec{u} = \nabla(\nabla \cdot \vec{u}) - \nabla \times (\nabla \times \vec{u})[/math]
[math]\nabla \cdot \vec{u}=\frac{\partial u_1}{\partial x}+\frac{\partial u_2}{\partial y} = \frac{\partial (y\cdot(1-y)\frac{p_1-p_2} {2μ})}{\partial x} +\frac{\partial 0}{\partial y} = 0[/math]
De este resultado se extrae que, al ser la divergencia nula, se cumple la condición de incompresibilidad. Esto implica que el fluido es incompresible (ni se comprime ni se expande) por lo que siempre ocupa el mismo volumen.
Para calcular el rotacional consideramos el campo en 3 dimensiones con z = 0:
[math]\nabla\times\vec u= \begin{vmatrix} \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{vmatrix}= \begin{vmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ \ (y-y^2)\frac{p_1-p_2}{2μ} & 0 & 0 \end{vmatrix}= -(1-2y)\frac{p_1-p_2}{2μ}\vec{k}[/math]
[math]\nabla\times(\nabla\times\vec u)= \begin{vmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ \ 0 & 0 & -(1-2y)\frac{p_1-p_2}{2μ} \end{vmatrix}= \frac{p_1-p_2}{μ}\vec{i}[/math]
Luego [math]\Delta\vec{u} = \frac{p_2-p_1}{μ}\vec{i}[/math]
Finalmente, sustituyendo en la ecuación:
Queda demostrado.
3 Campo de velocidades
3.1 Definición
El campo vectorial de velocidades objeto del estudio viene definido por la siguiente expresión:
[math]\vec{u}(x,y) = \frac{y(1-y)(p_1-p_2)}{2\mu}\vec{i}[/math]
Se observa fácilmente que este campo vectorial solamente guarda una relación de dependencia con la variable [math]y[/math], siendo independiente de los valores que le asignemos a [math]x[/math]. Esto se traduce en que la velocidad del fluido en un punto sólo depende de la distancia de este a las paredes laterales del canal. Además, la única componente del campo [math]\vec{u}(x,y)[/math] es la del vector [math]\vec{i}[/math], por lo que la velocidad será paralela a las paredes laterales independientemente del punto que escojamos.
3.2 Representación gráfica
Para representarlo gráficamente, suponemos:
[math]p_1=2 \qquad p_2=1 \qquad \mu=1[/math]
[math]\vec{u}(x,y) = \frac{y(1-y)(2-1)}{2}\vec{i}[/math]
Creamos a continuación el siguiente programa en MATLAB:
% Representación gráfica del campo de velocidades.
hold off
% Campo vectorial de velocidades:
[X,Y] = meshgrid([0:0.25:4],[0:0.1:1]);
ux=inline('y.*(1-y)./(2)','x','y');
uy=inline('0.*x','x','y');
U=ux(X,Y);
V=uy(X,Y);
quiver(X,Y,U,V)
% Ejes, región y rectas superior e inferior:
grid on
axis([0,4,-1,2])
axis equal
xlabel('X')
ylabel('Y')
hold on
x = [0,4];
y = [0,0];
plot(x,y,'-k','linewidth',2)
hold on
x = [0,4];
y = [1,1];
plot(x,y,'-k','linewidth',2)
3.3 Velocidad máxima
Los puntos correspondientes a la máxima velocidad del fluido, se determinarán igualando la primera derivada del campo de velocidades a cero:
[math]\vec{u}(x,y) = \frac{y(1-y)(p_1-p_2)}{2\mu}\vec{i}[/math]
[math]\frac{\partial }{\partial y} u(x,y) = \frac{(1-2y)(p_1-p_2)}{2\mu} = 0[/math]
[math]1-2y=0[/math]
[math]y=0.5[/math]
El resultado obtenido puede verificarse fácilmente observando la representación gráfica del campo de velocidades.
3.4 Rotacional del campo de velocidades
El rotacional del campo de velocidades es:
[math]\nabla\times\vec u= \begin{vmatrix} \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{vmatrix}= \begin{vmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ \ (y-y^2)\frac{p_1-p_2}{2μ} & 0 & 0 \end{vmatrix}= -(1-2y)\frac{p_1-p_2}{2μ}\vec{k}[/math]
Para la visualización del módulo del rotacional suponemos:
[math]p_1=2 \qquad p_2=1 \qquad \mu=1[/math]
[math]rot(x,y) =-(1-2y)/2[/math]
Creamos a continuación el siguiente programa en MATLAB:
% Representación gráfica del modulo del campo de rotacional.
hold off
% Modulo del rotacional:
[X,Y] = meshgrid([0:0.1:4], [0:0.1:1]);
Z = (-1+2*Y)/2;
surf(X,Y,Z)
% Ejes, región y rectas superior e inferior:
grid on
axis([0,4,-1,2])
axis equal
xlabel('X')
ylabel('Y')
colorbar
hold on
x = [0,4];
y = [0,0];
z = [-0.5,-0.5];
plot3(x,y,z,'-k','linewidth',2)
hold on
x = [0,4];
y = [1,1];
z = [0.5,0.5];
plot3(x,y,z,'-k','linewidth',2)
Como puede observarse en el gráfico, el módulo del rotacional únicamente depende de la variable [math]y[/math], siendo cte su valor para cualquier [math]x[/math]. Podemos ver que los puntos donde alcanza su máximo valor son aquellos donde la variable [math]y[/math] llega a los extremos de su intervalo [0,1], mientras que alcanza su valor mínimo en el punto intermedio del intervalo, en [math]y=0,5[/math]. ¿Son estos valores lógicos?
Si observamos el campo de velocidades del fluido vemos que alcanza sus valores máximos en los valores medios del intervalo, siendo la máxima velocidad del fluido la que se da en [math]y=0,5[/math]. Entendiendo por rotacional el campo que caracteriza el efecto de giro o rotación de un fluido, es comprensible que en los extremos del canal este efecto de giro *sea máximo en el sentido de desplazamiento el cual viene dado por el campo de velocidades, en nuestro caso el sentido positivo de la [math]x[/math]*.
Así pues el rotacional en [math]y=1[/math] tendrá valor máximo en el sentido de [math]\vec{k}[/math] positivo, lo que implica una rotación máxima de nuestro fluido en el sentido antihorario, mientras que en el otro extremo [math]y=0[/math] tendremos el mismo valor pero de signo negativo, induciendo rotación en sentido horario. *En ambos casos estos sentidos favorecen la circulación de nuestro fluido*. En valores de [math]y[/math] intermedios, donde la velocidad es máxima, el rotacional tiende 0 debido a que *la rotación del líquido es igual en ambos sentidos*.
3.5 Caudal del canal
El caudal es la cantidad de fluido que lleva el canal. Se obtiene observando la cantidad de fluido que atraviesa una determinada sección del mismo, por lo que al tratarse de un estudio en dos dimensiones sus unidades serán [math]\frac{m^2}{s}[/math]. Si queremos obtenerlo tendremos que calcular el flujo del campo de velocidades en una determinada superficie, en nuestro caso, al tratarse de un canal en dos dimensiones será el flujo que atraviese una curva. Se calculará resolviendo esta integral:
- [math] Caudal = \int_γ \vec {u(t)}\cdot \vec N \; ds [/math]
Siendo [math] \vec N \ [/math] un vector normal a la curva y [math]t[/math] el parámetro que la define. Si tomamos la curva de mínima distancia que cruza el canal, el caudal calculado será el total. El vector normal tendrá la misma dirección que el campo de velocidades y conservaremos la parametrización cartesiana para definir la curva en función de [math]y[/math], por lo que al caudal se obtendrá a partir de la integral::[math] \int_0^1 \frac{y(1-y)(p_1-p_2)}{2μ} \; dy = \frac{p_1-p_2}{2μ}\ \int_0^1 y(1-y) \; dy [/math] Resolveremos esta integral numéricamente.
El código matlab utilizado es el siguiente:
N=200;
a=0; b=1;
h=(b-a)/N;
M=0;
f=inline('y*(1-y)','y');
for i=1:N
y0=a+(i-1)*h;
y1=a+i*h;
M=M+(h/2)*(f(y0)+f(y1));
end
disp(M)
Obtenemos que el resultado de la integral es 0,1667 . Luego el caudal que circula por el canal será:
- [math] Caudal = \frac{(p_1-p_2)}{2μ}0,1667 [/math]
Suponiendo [math]p_1=2[/math], [math]p_2=1[/math] y [math]μ=1[/math] tendremos:
[math] Caudal = 0,083 \frac {m^2}{s} [/math]
3.5.1 Porcentaje de caudal central
Queremos conocer el porcentaje que refleje el caudal central con relación al total. Para ello calculamos el caudal central de la misma forma que antes solo que esta vez acotando la curva entre 0.25 y 0.75. Entonces si a=0.25 y b=0.75::[math] \frac{p_1-p_2}{2μ} \int_a^b y(1-y) \; dy [/math]
Resolvemos con el mismo programa matlab de antes solo que ahora cambiamos:
a=0.25; b=0.75;
Que resulta 0.1174. Por lo que el porcentaje del caudal que circula en la zona central es:
[math] Caudal_{central} = \frac{0.1174}{0.1667}100=70.4259 \% [/math]
Este resultado refleja que la velocidad es mayor en la zona central que en los bordes como ya observamos en la representación del campo.
4 Campo de presiones
4.1 Definición
El campo escalar de presiones objeto del estudio viene definido por la siguiente expresión:
[math]p(x,y) = p_1+(p_2-p_1)(x-1)[/math]
4.2 Representación gráfica
Para representarlo gráficamente, suponemos:
[math]p_1=2 \qquad p_2=1 \qquad \mu=1[/math]
[math]p(x,y) =3-x[/math]
Creamos a continuación el siguiente programa en MATLAB:
% Representación gráfica del campo de presiones.
hold off
%Campo escalar de presiones:
[X,Y] = meshgrid([0:0.1:4], [0:0.1:1]);
f=inline('3-x','x','y');
Z=f(X,Y);
surf(X,Y,Z)
%Ejes, región y rectas superior e inferior:
grid on
axis([0,4,-1,2])
axis equal
xlabel('X')
ylabel('Y')
colorbar
hold on
x = [0,4];
y = [0,0];
z = [3,-1];
plot3(x,y,z,'-k','linewidth',2)
hold on
x = [0,4];
y = [1,1];
z = [3,-1];
plot3(x,y,z,'-k','linewidth',2)
4.3 Presión media del fluido en el canal
Para hallar la presión media hemos empleado el Teorema del Valor Medio. La integral es fácilmente resoluble analíticamente, pero también la hemos resuelto numéricamente.
4.3.1 Resolución analítica
Sea [math]v[/math] el valor medio que pretendemos hallar, [math]D = [a,b]\times[c,d] = [0,4]\times[0,1][/math] el dominio que contiene el área que estamos estudiando del canal y A dicha área:
- [math] \frac{1}{A} \int_a^b \int_c^d f(x,y) \; dx \; dy = \frac{1}{4} \int_0^4 \int_0^1 p_1+(p_2-p_1)(x-1) \; dx \; dy = \frac{1}{4} \int_0^4 [ \int_0^1 p_1+(p_2-p_1)(x-1) \; dy \; ]dx = \frac{1}{4} \int_0^4 p_1+(p_2-p_1)(x-1) \; dx = [/math]
- [math] = \left. xp_1+(p_2-p_1)(\frac{x^2}{2}-x) \right |_0^4 \frac{1}{4} = \frac{1}{4} [(4p_1)+(p_2-p_1)(\frac{16}{2}-4)] = p_2[/math]
La presión disminuye linealmente en función de x desde 0 hasta 4, por lo que es lógico que el valor medio se encuentre justo en la mitad, en x = 2. Esto es coherente con lo que aparece en el gráfico.
4.3.2 Resolución numérica
Al aproximar la función por la Regla del Trapecio se obtiene un resultado exacto, al ser esta función lineal. El código utilizado es el siguiente:
N1=400; N2=100;
a = 0; b = 4; c= 0; d = 1;
h1= (b-a)/N1; h2=(d-c)/N2;
x=a:h1:b; y=c:h2:d;
[xx,yy] = meshgrid(x,y);
p1=2; p2=1;
f = p1 +(p2-p1)*(xx.-1);
w1 = ones(N1+1,1);
w1(1) = 1/2; w1(N1+1) = 1/2;
w2 = ones(N2+1,1);
w2(1) = 1/2; w2(N2+1) = 1/2;
result= (h1*h2*w2'*f*w1)/4
También obtenemos en esta ocasión [math]p_2[/math], teniendo en cuenta que [math]p_2 = 1[/math] en nuestro fluido.
5 Campo de temperaturas
5.1 Definición
La temperatura a lo largo del canal viene dada por el siguiente campo escalar:
[math]T(x,y)=e^{-(x-1)^2+y^2}[/math]
5.2 Representación gráfica
Para obtener su representacion nos ayudamos de Octave UPM:
% Representación gráfica del campo de temperaturas
hold off
% Creación del Mallado:
x=linspace(0,4,60);
y=linspace(0,1,15);
[Mx,My]= meshgrid(x,y);
% Campo escalar de temperaturas
Mz= exp(-(Mx-1).^2+My.^2);
% Creación de la superficie
surf(Mx,My,Mz)
% Curvas de nivel(color verde)
hold on
contour(Mx,My,Mz,25,'g','LineWidth',1)
% Ejes, región y rectas superior e inferior:
grid on
axis([0,4,-1,2])
axis equal
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')
colorbar
hold on
x = [0,4];
y = [0,0];
plot(x,y,'-k','linewidth',2)
hold on
x = [0,4];
y = [1,1];
plot(x,y,'-k','linewidth',2)
Comprobamos a continuación que el gradiente de temperaturas es ortogonal a las curvas de nivel de la temperatura.
Para ello escribimos el siguiente programa en MATLAB:
% Representación gráfica del campo de temperaturas.
hold off
% Mallado:
x=linspace(0,4,40);
y=linspace(0,1,10);
[X,Y]=meshgrid(x,y);
% Curvas de nivel y campo vectorial:
Z= exp(-(X-1).^2+Y.^2);
U=-2*(X-1).*exp(-(X-1).^2+Y.^2);
V=2.*Y.*exp(-(X-1).^2+Y.^2);
hold on
quiver(X,Y,U,V)
contour(X,Y,Z,25);
% Ejes, región y rectas superior e inferior:
grid on
axis([0,4,-1,2])
axis equal
xlabel('X')
ylabel('Y')
hold on
x = [0,4];
y = [0,0];
plot(x,y,'-k','linewidth',2)
hold on
x = [0,4];
y = [1,1];
plot(x,y,'-k','linewidth',2)
Ampliando al recinto [1.25; 1.75] x [0.25; 0.75] se puede apreciar mejor la ortogonalidad:
6 Líneas de corriente
6.1 Definición
Las líneas de corriente (líneas tangentes a [math] \vec{u}(x,y) [/math] en cada punto) vienen definidas por las isolíneas del potencial [math]\psi[/math] del campo [math]\vec{v}(x,y)[/math], siendo el campo [math] \vec{v}(x,y) = \vec{k}\times \vec{u} [/math] perpendicular al campo [math] \vec{u}(x,y) [/math]
[math]\vec{v}(x,y) = \vec{k} \times \vec{u}(x,y) = \frac{y(1-y)(p_1-p_2)}{2\mu}\vec{j}= \frac{1}{2}y(1-y)\vec{j}= \frac{1}{2}(y-y^2)\vec{j}[/math]
[math]\psi_v=\frac{1}{2}(\frac{y^2}{2}-\frac{y^3}{3})[/math]
6.2 Representación gráfica
Creamos a continuación el siguiente programa en MATLAB:
% Representación gráfica de las isolíneas del potencial del campo
% v = k^u, correspondientes a las lineas de corriente del campo u.
hold off
% Potencial 'Ψ' correspondiente al campo v:
[X,Y] = meshgrid([0:0.5:4], [0:0.05:1]);
p=inline('(0.5).*(((y.^2)./2)-((y.^3)./3))','x','y');
P=p(X,Y);
contour(X,Y,P,15,'-r')
% Ejes, región y rectas superior e inferior:
grid on
axis([0,4,-1,2])
axis equal
xlabel('X')
ylabel('Y')
hold on
x = [0,4];
y = [0,0];
plot(x,y,'-k','linewidth',2)
hold on
x = [0,4];
y = [1,1];
plot(x,y,'-k','linewidth',2)
Tal como se puede observar en la representación gráfica, las isolíneas del potencial del campo v(x,y) → ‘Ψv’, corresponden a las lineas de corriente del fluido.
