Visualización de campos vectoriales y escalares en fluidos (grupo C6)
Contenido
1 Introducción
Consideraremos el flujo de un fluido incompresible a través de un canal con paredes rectas que representaremos (en coordenadas cartesianas) sobre el siguiente mallado.
x0=0;xN=4;dx=(xN-x0)/40;
y0=0;yM=1;dy=(yM-y0)/10;
x=x0:dx:xN;
y=y0:dy:yM;
[xx,yy]=meshgrid(x,y);
figure(100)
mesh(xx,yy,0.*xx)
axis([0,4,-1,2])
view(2)
Al tratarse de un fluido incompresible, es decir, el flujo ni se expande ni se contrae, luego cumplirá que la divergencia será cero y satisfará la ecuación de Navier-Stokes:
La velocidad de las partículas de nuestro fluido viene dada por [math]\vec u(x,y)= y(1-y)(p_1-p_2)/(2\mu) \vec i[/math] y su presión [math]p(x,y)=p_1 + (p_2-p_1)(x-1)[/math], donde [math]p_1[/math] es la presión en los puntos [math]x=1[/math], [math]p_2[/math] la presión en los puntos x=2 y [math]\mu [/math] el coeficiente de viscosidad del fluido. Comprobaremos que ([math]\vec u[/math],p) satisface las condiciones anteriores.
si [math]\nabla \cdot \vec{u}=0 [/math] , entonces [math]\vec{u} \cdot \nabla \vec{u} = 0 [/math]
[math]\nabla p=μ \Delta\vec{u} [/math]
[math]\nabla p=dp_i/dx_i =p_2-p_1[/math]
[math]\Delta\vec{u} = \displaystyle\frac{d^2u}{dx^2}+\displaystyle\frac{d^2u}{dy^2}+\displaystyle\frac{d^2u}{dz^2} = \displaystyle\frac{d^2u}{dy^2}[/math]
[math]\displaystyle\frac{d^2u}{dy^2} = \displaystyle\frac{p_2-p_1}{\mu}[/math]
Como [math]\nabla p=μ \Delta\vec{u} [/math], entonces [math]p_2-p_1 = \mu( \displaystyle\frac{p_2-p_1}{\mu})[/math] quedando demostradas las condiciones iniciales de nuestro fluido incompresible.
Ahora comprobaremos que la velocidad del fluido [math]\vec u(x,y)= y(1-y)(p_1-p_2)/(2\mu) \vec i[/math] en las paredes (y = 0, y =1) es nula, debido a la fricción que se genera entre las paredes y las partículas del fluido en esa zona.
En y = 0
[math]\displaystyle\frac{0(1-0)(p_1-p_2)}{2\mu} = \displaystyle\frac{0}{2\mu} = 0[/math]
En y = 1
[math]\displaystyle\frac{1(1-1)(p_1-p_2)}{2\mu} = \displaystyle\frac{0}{2\mu} = 0[/math]
2 Campo de velocidades y presiones
Suponiendo que [math]p_1=2[/math], [math]p_2=1[/math] y [math]\mu=1[/math] el campo de presiones y el de velocidades es
[math]\vec u(x,y)=\displaystyle\frac{(y-y^2)}{2}\vec i[/math] [math] p(x,y) = 3-x[/math]
close all
hold on
q=linspace(0,4,100);
p=linspace(0,1,10);
[Q,P]=meshgrid(q,p);
Z=3-Q;
pcolor(Q,P,Z);
axis([0,4,-1,2])
x=linspace(0,4,10);
y=linspace(0,1,10);
[X,Y]=meshgrid(x,y);
U=(Y.*(1-Y))/2;
V=U.*0;
quiver(X,Y,U,V,'k','Linewidth',1.5)
axis([0,4,-1,2])
hold off
Como se ha comprobado analíticamente y ahora podemos comprobar gráficamente, la velocidad del fluido es máxima en el centro y nula en las paredes. Además, también podemos comprobar que la presión es constante en la abscisa y disminuye conforme nos alejamos del origen.
3 Líneas de corriente del campo [math]\vec{u} [/math]
Dibujaremos ahora las líneas de corriente del campo [math]\vec u[/math], es decir, las líneas que son tangentes a [math]\vec u[/math] en cada punto. Para ello calculamos el campo [math]\vec v[/math] que en cada punto es ortogonal a [math]\vec u[/math] para ello tomaremos [math]\vec v=\vec k \times \vec u[/math].
[math]\vec v=\vec k \times \vec u = (\displaystyle\frac{p_1-p_2}{2\mu})(y-y^2)\vec j[/math]
[math]\nabla \times \vec v = \frac{y(1-y)}{2}\vec{j} = (\nabla·\vec u)\vec k = \vec 0[/math]
Como el campo [math]\vec{u} [/math] es adivergente, es decir, su divergencia es 0, conlleva a que el campo [math]\vec{v} [/math], ortogonal a [math]\vec{u} [/math], sea irrotacional, ya que al calcular el rotacional es igual a [math](\nabla·\vec u)\vec k [/math].
Vamos a calcular el potencial escalar [math]\psi[/math], también conocido como función de corriente de [math]\vec{u} [/math]
[math]\vec{v} = \nabla\psi = (\displaystyle\frac{d\psi }{dx})\vec{i} + (\displaystyle\frac{d\psi }{dy})\vec{j} [/math] por lo que
[math]\psi(x,y) = \displaystyle\frac{p_1-p_2 }{2\mu}(\displaystyle\frac{y^2}{2}-\displaystyle\frac{y^3}{3})+cte[/math]
[math]\frac{ \partial \psi}{\partial y}=\frac{y(1-y)}{2} [/math] ; [math]\psi =\int\frac{y(1-y)}{2}dy =\frac{1}{2}(\frac{y^2}{2}-\frac{y^3}{3})+ f(x) [/math]
[math]\frac{ \partial \psi}{\partial x}=0 [/math] ; [math]\frac{\partial}{\partial x}(\frac{1}{2}(\frac{y^2}{2}-\frac{y^3}{3})) + f'(x)=0 \implies f'(x)=0 \implies f(x)=cte[/math]
hold on
x=linspace(0,4,10);
y=linspace(0,1,10);
[X,Y]=meshgrid(x,y);
U=(Y.*(1-Y))/2;
V=U.*0;
quiver(X,Y,U,V,'k','Linewidth',1)
axis([0,4,-1,2])
q=linspace(0,4,10);
p=linspace(0,1,10);
[xx,yy]=meshgrid(q,p);
psi= (1/2).* ((yy.^2/2)-(yy.^3/3));
contour(xx,yy,psi)
axis([0,4,-1,2])
hold off
Ahora vamos a calcular los puntos en los que la velocidad del fluido es máxima, para ello derivamos el campo de velocidades e igualamos a cero. Sabiendo que [math] \vec u(x,y) = \displaystyle\frac{1}{2}(y-y^2)\vec i = u(y)\vec i[/math] entonces
[math]\displaystyle\frac{du}{dy} = 1-2y = 0[/math] por lo que la velocidad máxima se alcanzará en [math] y = \displaystyle\frac{1}{2}[/math]
[math] y = \displaystyle\frac{1}{2}[/math] corresponde a los puntos equidistantes de las paredes del canal. Esto se debe a que en las paredes existe rozamiento, que acaba anulando la velocidad de las partículas del fluido en contacto con las paredes, como hemos visto anteriormente.
4 Rotacional [math] \vec u [/math]
Primero vamos a calcular el rotacional de [math]\vec{u}[/math], que es el operador vectorial que muestra la tendencia de un campo vectorial a rotar alrededor de un punto.
[math]\nabla\times\vec{u}= \begin{pmatrix} \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{pmatrix}= \begin{pmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ \frac{y(1-y)(p_1-p_2)}{2\mu}& 0 & 0 \end{pmatrix}=-\frac{\partial}{\partial y}(\frac{y(1-y)(p_1-p_2)}{2\mu})\vec{k} =-\frac{ (1-2y)(p_1-p_2)}{2\mu}\vec{k} [/math]
Ahora particularizamos a los valores de [math] p_1=2[/math], [math] p_2=1[/math] y [math] \mu=1 [/math],
Hallamos su módulo:
El siguiente paso es ver cuando el rotacional del campo alcanza su máximo y mínimo. Para lo cual calculamos la derivada del modulo y la igualamos a 0:
como vemos, la derivada es una constante, luego no se alcanzan ni máximos ni mínimos. El máximo y mínimo lo obtenemos en los extremos del intervalo, en las paredes de la tubería, ya que al ser en estos puntos la velocidad nula, la partícula tenderá a su máxima rotación alrededor de ese punto en toda la zona. De la misma manera, al tener en el centro velocidad máxima, indica que el rotacional es cero, es decir, la partícula no tenderá a girar sobre los puntos situados en y=1/2.
[math]y=0[/math]
[math]|\nabla\times\vec{u}|=|\frac{ (1-2y)}{2}|=|\frac{1}{2} |=\frac{1}{2}[/math]
[math]y=1[/math]
[math]|\nabla\times\vec{u}|=|\frac{ (1-2y)}{2}|=|- \frac{1}{2}| =\frac{1}{2}[/math]
x0=0;xN=4;dx=(xN-x0)/40;
y0=0;yM=1;dy=(yM-y0)/10;
x=x0:dx:xN;
y=y0:dy:yM;
[xx,yy]=meshgrid(x,y);
p1=2;p2=1;mu=1;
pp=p1+(p2-p1)*(xx-1);
ff=((p1-p2)/(2*mu))*((1/2)*yy.^2-(1/3)*yy.^3);
Ux=((p1-p2)/(2*mu))*(yy-yy.^2);
Uy=0.*xx;
Vx=0.*xx;
Vy=((p1-p2)/(2*mu))*(yy-yy.^2);
rot=(.5-yy);
rot1=abs(.5-yy);
hold on
contour(xx,yy,rot,12);
axis([0,4,-1,2])
hold off
hold on
surf(xx,yy,rot);
axis([0,4,-1,2])
hold off
figure(660)
subplot(1,2,1)
surf(xx,yy,rot);
subplot(1,2,2)
surf(xx,yy,rot1);
5 Campo de temperaturas
Supongamos que la temperatura del fluido viene dada por el campo [math]T(x,y)=(x-1)^2-y^2[/math]. En la siguiente gráfica se representa el campo de temperaturas y su gradiente donde podremos comprobar que este es perpendicular a las curvas de nivel de la temperatura. Esto se debe a que el gradiente es la dirección de máximo crecimiento de la temperatura en cada punto.
x=x0:dx:xN;
y=y0:dy:yM;
[xx,yy]=meshgrid(x,y);
TT=(xx-1).^2-yy.^2;
gradTx=2*(xx-1);
gradTy=-2*yy;
figure(700)
hold on
contour(xx,yy,TT,40)
quiver(xx,yy,gradTx,gradTy)
axis([0,4,-1,2])
view(2)
hold off
