Visualización de campos en fluidos (Grupo C12)
Contenido
Las ecuaciones de Navier-Stokes son un conjunto de ecuaciones que describen el movimiento de fluidos. Con ellas se pueden analizar corrientes oceánicas, flujo alrededor de aviones o vehículos o cualquier otro fluido newtoniano presente en la atmósfera de la Tierra. Estas ecuaciones son en derivadas parciales no lineales, y por ello no se tiene solución para las mismas. La mayoría de las veces no se pueden usar analíticamente, salvo ciertos tipos de flujo y situaciones muy concretas, y se tiene que llevar a cabo un análisis numérico para obtener una solución aproximada.
Las ecuaciones de Navier-Stokes para un fluido son la estacionaria (1) y la condición de incompresibilidad (2), que se comprueban a continuación:
Siendo $\vec{u}$ la velocidad del fluido y p la presión del mismo.
Sea un fluido con velocidad [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] y presión [math]\vec{p}(x,y)=p_1+(p_2-p_1)(x-1)[/math] donde $p_1$ y $p_2$ son la presión en x=1 y x=2 respectivamente
1.1 Ecuación estacionaria
Se comprueba si el fluido satisface la ecuación de Navier-Stokes estacionaria
Esta expresión en lenguaje indicial será:
Al trabajar en coordenadas cartesianas las componentes contravariantes y covariantes coinciden. Ahora:
Operando
[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]
[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]
[math]\begin{pmatrix} (\Delta\vec{u})_1 \\ (\Delta\vec{u})_2 \end{pmatrix} = \begin{pmatrix} \frac{\partial^2 u_1}{\partial x^2}+\frac{\partial^2 u_1}{\partial y^2} \\ \frac{\partial^2 u_2}{\partial x^2}+\frac{\partial^2 u_2}{\partial y^2} \end{pmatrix} = \begin{pmatrix} \frac{p_2-p_1}{μ} \\ 0 \end{pmatrix}[/math]
Se comprueba si se cumple la igualdad
Queda comprobado.
1.2 Condición de incompresibilidad
La comprobación de las ecuaciones de Navier-Stokes pasa ahora por verificar que se trata de un fluido incompresible, es decir, que ni se expande ni se contrae:
[math]\nabla \cdot \vec{u}=0 [/math]
La divergencia de $\vec{u}$ se define como:
[math]\nabla \cdot \vec{u}=\frac{\partial u_1}{\partial x}+\frac{\partial u_2}{\partial y} [/math]
En el caso estudiado, [math]u_2(x,y)=0[/math] y, por tanto:
[math]\nabla \cdot \vec{u}=\frac{\partial u_1}{\partial x}=\frac{\partial}{\partial x}(\frac{y(1-y)(p_1-p_2)}{2μ})=0[/math]
2 Campo de presiones y campo de velocidades
Para el estudio de los campos de presiones y velocidades, consideramos el flujo del fluido a través de un canal con paredes rectas de dimensiones 1mx4m.
Trabajaremos en el plano con coordenadas cartesianas. Empezamos definiendo un mallado de ejes [0,4]x[-1,2] y dentro de ellos el fluido que ocupa el rectángulo [0,4]x[0,1], simulando el canal. El codigo matlab para realizarlo es el siguiente
x=0:0.1:4; %Mallado en el intervalo [0,4]x[0,1] con longitud de paso de 0,1
y=0:0.1:1;
[xx,yy]=meshgrid(x,y); %Hacemos el mallado
mesh(xx,yy,xx*0)
axis([0,4,-1,2]) %Definimos los ejes
2.1 Campo de presiones
Veamos el campo de presiones del fluido:
Tomando [math] p_{1}=2[/math], [math] p_{2}=1[/math] y [math] \mu=1[/math] los campos de presiones y velocidades que resultan:
%Dibujo del campo
x=0:0.1:4; %Mallado en el intervalo [0,4]x[0,1] con longitud de paso de 0,1
y=0:0.1:1;
[xx,yy]=meshgrid(x,y);
%Sólido antes de desplazarse
figure(1);
subplot(1,2,1)
p= 3 - xx;
surf(xx,yy,p)
%axis([0,4,-1,2])
figure(2)
subplot(1,2,1)
p=3-xx;
surf(xx,yy,p)
axis([0,4,-1,2])
El campo de presiones viene representado por una ecuación, la cual sólo depende de la variable x, lo que nos conduce a afirmar que para todo valor de y arbitrario, la presión depende de x exclusivamente. Esto explica y da sentido a nuestras gráficas ya que la zona de color rojo nos indica altas presiones al ser los valores de x próximos a cero, y alcanzar la función máximos valores para cualquier valor de y. A medida que aumentamos los valores de x la presión disminuye y alcanza ese tono azul.
Para obtener la gráfica en 2D simplemente proyectamos.
2.2 Campo de velocidades
%visualización del campo de velocidades
x=0:0.5:4; %Mallado en x
y=0:0.1:1; %Mallado en y
[xx,yy]=meshgrid (x,y);
ux=1/2.*(yy-yy.^2);
uy=0.*xx
figure(1)
quiver(xx,yy,ux,uy)
2.2.1 Puntos de velocidad nula y velocidad máxima
A continuación, se investiga en qué puntos la velocidad es nula y en cuáles es máxima. En primer lugar, se iguala la expresión de la velocidad al vector $\vec{0}$, para averiguar dónde es nula
[math]\vec{u}=y(1-y)\frac{p_1-p_2}{2μ}\vec{i}=\vec{0}[/math]
[math] y(1-y)=0 \left \{ \begin{matrix} y=0 \\ y=1 \end{matrix}\right. [/math]
Para hallar los puntos de velocidad máxima, se obtienen los extremos relativos de la función [math]g(y)=y-y^2[/math], que serán los mismos que los de la función dada.
[math]g'(y)=0 \longrightarrow 1-2y=0 \longrightarrow y=0.5[/math]
[math]g''(y)=-2\lt0[/math]
La función presenta un máximo en y=0.5
Corroboramos lo visto analíticamente mediante la siguiente gráfica:
%Dibujo del campo
x=0:0.1:4; %Mallado en el intervalo [0,4]x[0,1] con longitud de paso de 0,1
y=0:0.1:1;
[xx,yy]=meshgrid(x,y);
%Sólido antes de desplazarse
figure(1);
subplot(1,2,1)
mesh(xx,yy,0*xx)
axis([0,4,-1,2])
%Sólido después de desplazarse
subplot(1,2,2)
ux=(1/2)*(yy.*(1-yy));
uy=0*xx;
mesh(xx+ux,yy+uy,0*xx)
axis([0,4,-1,2])
En la primera gráfica representamos el sólido antes de desplazarse en dicho mallado de puntos, y proyectado según los ejes.
Para representar el desplazamiento longitudinal, usamos el campo de velocidades del fluido, el cual depende de la variable y (altura del punto). Sabiendo que s=vt y para tiempo igual a 1 segundo, s=v, o lo que es lo mismo, el valor del desplazamiento coincide con el de la velocidad. Por lo tanto, sumamos al punto su velocidad y vemos el desplazamiento que experimenta. De esta forma, en las paredes del canal, no hay desplazamiento por tanto la velocidad es nula y en el centro, al tener el máximo valor de y, es máxima. Esto es lo que queda representado en la segunda gráfica.
La tercera gráfica muestra todo lo anterior pero en una sección x.
2.2.2 Líneas de corriente
Vamos a dibujar ahora las líneas de corriente del campo de velocidades $\vec{u}$, es decir, las líneas de corriente son tangentes a $\vec{u}$en cada punto. De forma general, siendo $\vec{u}$ un campo cualquiera que tenga Ψ como potencial escalar, las superficies equipotenciales (superficies con el mismo potencial) se caracterizan por ser ortogonales al campo en cada punto. En el caso tratado, al trabajar en $\mathbb{R^2}$ , las superficies equipotenciales pasarán a ser líneas equipotenciales y se definirán como la intersección entre la superficie equipotencial y el plano en el que se trabaja (en nuestro caso el plano XOY) Para definir las líneas de corriente $\vec{u}$ (tangentes al campo $\vec{u}$ en cada punto) vamos a definir un campo $\vec{v}$ ortogonal al mismo en cada punto y definir éstas como las líneas equipotenciales de $\vec{v}$.
[math] \vec{v}=\vec{k}\times \vec{u} = \displaystyle\frac{y(1-y)(p_1-p_2)}{n}\vec{j}[/math]
Considerando
[math]\vec{v}: [0,4]x[0,1]\subset\mathbb{R^2}\longrightarrow\mathbb{R^2}[/math]
[math](x,y)\longrightarrow\vec{v}(x,y)=\displaystyle\frac{y(1-y)(p_1-p_2)}{n}\vec{j}[/math]
Existe
[math]\exists[0,4]x[0,1]: \mathbb{R^2}\longrightarrow\mathbb{R}[/math]
tal que
[math]\nablaΨ\longrightarrow\vec{v}(x,y)[/math]
[math] \left . \begin{matrix} \frac{\partial Ψ}{\partial x}=u_1(x,y) \\ \frac{\partial Ψ}{\partial y}=u_2(x,y) \end{matrix} \right \}[/math]
Del sistema de ecuaciones diferenciales se deduce el potencial Ψ
Su representación gráfica en matlab es:
x=0:0.25:4;
y=0:0.2:1;
[xx,yy]=meshgrid(x,y); % Mallado en [0,4]x[0,1] con equidistancia 0.1
p=(1/2).*((yy.^2/2)-(yy.^3/3)); % Potencial escalar
contour(xx,yy,p) %Líneas de nivel
axis([0,4,-1,2]) % Ejes en [0,4]x[-1,2]
2.2.3 Rotacional
[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]
¿Es razonable?
Entendamos el rotacional como un vector perpendicular al plano de curvatura del campo que caracteriza el efecto de giro o rotación del fluido. Para entender esto podemos verlo con el ejemplo práctico de "la paleta":
Para ver la rotación en un punto, cogemos una sección a la que pertenezca el punto y dibujamos su campo de velocidades. Como éste es constante en x, es lógico coger una sección cualquiera "x". (ver figura adjunta). Nos imaginamos al lado del campo de velocidades una "paleta" que tiene como eje de giro el punto en cuestión. El campo de velocidades "empujará" la paleta de una forma u otra en función de dónde se encuentre dicho eje. Veamos algunos ejemplos: (EJEMPLO EN MEDIO, EN LA MITAD SUPERIOR Y EN EXTREMO)
Si el eje está en medio, el campo por arriba y por abajo son iguales, pero empujan en sentido contrario por tanto, el rotacional será nulo.
Si el eje está en el medio de la parte superior, girará en sentido antihorario porque el empuje del campo por abajo es mayor que el de arriba.
Por último, si el eje está en el extremo superior, todo el campo empujará la paleta en sentido antihorario.
En la parte inferior sería totalmente análogo al ser simétrico.
%Valor absoluto del rotacional
y=0:0.1:1;
x=abs(1/2*(1-2*y));
plot(x,y)
3 Campo de temperatura
Supongamos que la temperatura del fluido viene dado por el campo: [math]T(x,y)=(x-1)^2-y^2 [/math]
Veamos su representación gráfica:
%Dibujo del campo
x=0:0.1:4; %Mallado en el intervalo [0,4]x[0,1] con longitud de paso de 0,1
y=0:0.1:1;
[xx,yy]=meshgrid(x,y);
%Sólido antes de desplazarse
figure(1);
subplot(1,2,1)
p=((xx-1).^2)-yy.^2;
surf(xx,yy,p)
%axis([0,4,-1,2])
figure(2)
subplot(1,2,1)
p=((xx-1).^2)-yy.^2
surf(xx,yy,p)
axis([0,4,-1,2])
Sea $\vec{w}$ el gradiente de la temperatura. Comprobemos gráficamente que son ortogonales.
f=inline('((x-1).^2)-y.^2','x','y'); %Campo escalar de temperatura
fx=inline('2.*(x-1)','x','y'); %derivadas parciales del campo
fy=inline('2.*y','x','y');
x=0:0.1:4;
y=0:0.1:1;
[X,Y]=meshgrid(x,y); %hacemos el mallado
Z=f(X,Y);
contour(X,Y,Z,20,'k'); %Curvas de nivel
hold on
xx=0:0.25:4;
yy=0:0.1:1;
[XX,YY]=meshgrid(xx,yy) %Representación de las flechas
U=fx(XX,YY);
V=fy(XX,YY);
quiver(XX,YY,U,V)
axis equal
OBSERVACIÓN: Esta propiedad, aplicada al campo de velocidades, es la utilizada para definir las líneas de corriente del mismo.
 Este artículo está en versión beta. El autor de este artículo no lo ha terminado todavía, por favor no lo edites hasta que elimine este mensaje.
|
.
