Diferencia entre revisiones de «Trabajo De Campos Grupo 7»

Revisión actual del 19:19 9 dic 2019

1 . Introducción

Trabajaremos con campos escalares y vectoriales en fluidos. En nuestro caso, vamos a considerar el flujo de un fluido incompresible a través de un canal con paredes rectas.

2 . Mallado que representa los puntos ocupados por un fluido

En primer lugar comenzaremos dibujando un mallado que represente los puntos interiores del rectángulo [0, 8] × [−1, 1] ocupado por un fluido. Fijaremos los ejes en la región [0, 4] × [−2, 2].

Mallado Inicial

%Definimos los intervalos
x=0:0.1:8
y=-1:0.1:1
[xx,yy]=meshgrid(x,y)
figure(1)
mesh(xx,yy,0*xx)
axis([0,4,-2,2])
view(2)

3 . Ecuación de Navier-Stokes

La velocidad de las partículas del fluido viene dada por:

y su presión por:

. A continuación verificaremos la ecuación de Navier-Stokes:

[math] \vec u \cdot \nabla \vec u + \nabla p = \mu Δ \vec u [/math]

Primero calculamos [math] \nabla \vec u [/math]:

[math] \nabla \vec u = \frac{ \partial u_{i}}{ \partial x^{j} } \vec g^{j} = 0\vec i + \frac{ \partial u_{1}}{ \partial y }\vec j = \frac{(-2y)(p_{1}-p_{2})}{(2\mu)}\vec j [/math]

[math] \vec u \cdot \nabla \vec u = \frac{(-2y)(p_{1}-p_{2})}{(2\mu)} \vec j \cdot \frac{(y+1)(1-y)(p_{1}-p_{2})}{(2\mu)} \vec i = 0 [/math]

El gradiente del campo es ortogonal al propio campo, por lo que su producto escalar es 0. Segundo calculamos [math] \nabla p [/math]:

[math] \nabla p = \frac{ \partial p}{ \partial x^{j} } \vec g^{j} = \frac{ \partial p }{ \partial x } \vec i + \frac{ \partial p }{ \partial y} \vec j = (p_{2}-p_{1}) \vec i [/math]

Por último, calculamos [math] Δ \vec u [/math]

[math] Δ \vec u = \nabla \cdot \nabla \vec u = \frac{ \partial^2 u_{i}}{ \partial^2 x^{j} } \vec g^{i} = \frac{ \partial^2 u_{1}}{ \partial y^2 } \vec i + 0 \vec j = (p_{2}-p_{1}) \vec i [/math]

Queda verificada la ecuación de Navier-Stokes, por lo que el fluido sigue un régimen estacionario.

Además se verifica la condición de incomprensibilidad .

Comprobamos que también la velocidad del fluido en las paredes es nula

4 . Campo de Presiones y Campo de Velocidades

Suponemos que el campo de presiones toma los siguientes valores:

[math] p_{1} = 2 [/math], [math] p_{2} = 1 [/math] y [math] \mu = 1 [/math]:

El campo de presiones viene dado por la expresión . Por otro lado, la expresión de las velocidades de las partículas viene dada por:

Sustituyendo los valores que nos dan, tenemos que:

El código para obtener la gráfica del campo de presiones es el siguiente:

Campo de Presiones

x=0:0.1:8;
y=-1:0.1:1;
[xx,yy]=meshgrid(x,y);
figure(1)
p=2+(1-2)*(xx-1);
surf(xx,yy,p)
axis([0,4,-2,2])
view(2)
colorbar
title('CAMPO DE PRESIONES')

Por último, para obtener el campo de velocidades y su gráfica, utilizamos el siguiente código:

Campo de Velocidades

x=0:0.1:8;
y=-1:0.1:1;
[xx,yy]=meshgrid(x,y);
ux=((yy+1).*(1-yy).*(2-1))./(2);
uy=0.*yy;
hold on
quiver(xx,yy,ux,uy)
axis([0,4,-2,2])
hold off
view(2)
title('CAMPO DE VELOCIDADES')

5 . Lineas de Corriente del Campo ū y su Función de Corriente

Para obtener las líneas de corriente, calculamos el siguiente producto vectorial:

x=0:0.1:8;
y=-1:0.1:1;
[xx,yy]=meshgrid(x,y);
figure(1)
ux=((yy+1).*(1-yy).*(2-1))./(2);
uy=0.*yy;
vy=((yy+1).*(1-yy).*(2-1))./(2);
vx=0.*yy;
hold on
quiver(xx,yy,ux,uy)
quiver(xx,yy,vx,vy)
axis([0,4,-2,2])
view(2)
hold off

6 . Puntos con Velocidad Máxima del Fluido

Los puntos con velocidad máxima del fluido serán aquellos cuya primera derivada de su módulo sea cero, y la segunda derivada en el que sea el punto de mayor velocidad, sea inferior a cero.

[math] \vec u (x,y) = \frac{(y+1)(1-y)(p_{1}-p_{2})}{(2\mu)}\vec i [/math]

[math] | \vec u (x,y) | = \frac{(y+1)(1-y)(p_{1}-p_{2})}{(2\mu)} [/math]

[math] \frac{d| \vec u (x,y) |}{dy} = \frac{(-2y)(p_{1}-p_{2})}{(2\mu)} [/math]

[math] \frac{d| \vec u (x,y) |}{dy} = 0 \implies \frac{(-2y)(p_{1}-p_{2})}{(2\mu)} = 0 \implies y = 0 [/math]

[math] \frac{d^2| \vec u (x,y) |}{dy^2} = \frac{-2(p_{1}-p_{2})}{(2\mu)} [/math]

Suponemos el máximo o mínimo en función de los valores [math] p_{1} = 2 [/math], [math] p_{2} = 1 [/math] y [math] \mu = 1 [/math]:

[math] \frac{d^2| \vec u (x,y) |}{dy^2} = \frac{-2 \cdot(2-1)}{(2\cdot 1)} = -1 [/math]

Sería menor de cero, por lo que la velocidad sería máxima en todos los puntos cuya ordenada sea [math] y = 0 [/math].

7 . Rotacional de ū

Vamos a calcular el rotacional del campo [math] \vec u [/math] y después su módulo: [math]\nabla\times\vec{u}= det\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}=det\begin{pmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ \frac{(y+1)(1-y)(p_1-p_2)}{2\mu}& 0 & 0 \end{pmatrix}=-\frac{\partial}{\partial y}(\frac{(y+1)(1-y)(p_1-p_2)}{2\mu})\vec{k} =\frac{ -(2y)(p_1-p_2)}{2\mu}\vec{k} [/math] [math]|\nabla\times\vec{u}| = \frac{(-2y)(p_1-p_2)}{2\mu}[/math] [math]\frac{d|\nabla\times\vec{u}|}{dy} = \frac{-2(p_1-p_2)}{2\mu} = \frac{-(p_1-p_2)}{\mu} \implies\frac{-(p_1-p_2)}{\mu}= 0 [/math]

No existen máximos ni mínimos relativos. Vamos a ver en las paredes si son máximos o mínimos absolutos en el intervalo [-1,1]:

[math] y = -1 \implies |\nabla\times\vec{u}| = |\frac{(-2\cdot{(-1)})(p_1-p_2)}{2\mu}| = |\frac{p_1-p_2}{\mu}| [/math]

[math] y = 1 \implies |\nabla\times\vec{u}| = |\frac{(-2\cdot{1})(p_1-p_2)}{2\mu}| = |\frac{-(p_1-p_2)}{\mu}| [/math]

Para [math] p_{1} = 2 [/math], [math] p_{2} = 1 [/math] y [math] \mu = 1 [/math]

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

[math] y = 1 \implies |\nabla\times\vec{u}| = |\frac{-(2-1)}{1}| = |-{1}|=1 [/math]

Los puntos de mayor rotacional se encuentran en las paredes del canal. Representamos el módulo del rotacional del campo [math] \vec u [/math], que sería la función [math]|\nabla\times\vec{u}| = {y} [/math]. Programando en Matlab:

Rotacional del Campo de Velocidades

x=0:0.1:8;
y=-1:0.1:1;
[xx,yy]=meshgrid(x,y);
mod_rot=abs(yy);
surf(xx,yy,mod_rot)
colorbar
view(2)
axis([0,4,-2,2])

8 . Temperatura del Fluido

La temperatura del fluido viene dada por la siguiente ecuación:

Como el Campo de Temperaturas viene en coordenadas polares, hemos de cambiarlo a coordenadas cartesianas. Definiremos la variable rho como: Por otro lado, teta será igual a Le sumamos a x un valor que tienda a -∞ para asegurarnos que el denominador de la fracción no sea igual a cero evitando que el programa nos de un error.

Representaremos el campo de temperaturas utilizando Matlab, con el siguiente código:

Temperatura del Fluido

x=linspace(0,8,30);
y=linspace(-1,1,30);
rho=sqrt(x.^2+y.^2)
teta=atan(y./[x+0.0000001])
[rho,teta]=meshgrid(x,y)
figure(1)
T=1+[sin(teta).^2].*[exp(-[rho-0.5].^2)]
surf(rho,teta,T)
hold on
plot(rho,teta,'k','linewidth',1)
axis([0,4,-2,2])
view(2)
contour(rho,teta,T,100)

Las curvas de nivel de la temperatura, se obtendrán gráficamente con el siguiente código:

Curvas de nivel de la temperatura

x=linspace(0,8,30);
y=linspace(-1,1,30);
rho=sqrt(x.^2+y.^2)
teta=atan(y./[x+0.0000001])
[rho,teta]=meshgrid(x,y)
T=1+[sin(teta).^2].*[exp(-[rho-0.5].^2)]
Tx=-2.*[sin(teta).^2].*[rho-0.5].*[exp(-[rho-0.5].^2)]
Ty=sin(2*teta).*[exp(-[rho-0.5].^2)]
figure(1)
contour(rho,teta,T,100)

9 . Gradiente de la Temperatura

A continuación hallaremos el gradiente de la temperatura utilizando la siguiente fórmula.

Por lo que haciendo las derivadas parciales obtenemos que el gradiente es:

Para representarlo en un gráfico utilizaremos el siguiente código Matlab:

Gradiente de la temperatura

x=linspace(0,8,30);
y=linspace(-1,1,30);
rho=sqrt(x.^2+y.^2)
teta=atan(y./[x+0.0000001])
[rho,teta]=meshgrid(x,y)
T=1+[sin(teta).^2].*[exp(-[rho-0.5].^2)]
%Primer componente del gradiente:
Tx=-2.*[sin(teta).^2].*[rho-0.5].*[exp(-[rho-0.5].^2)]
%Segundo componente del gradiente:
Ty=sin(2*teta).*[exp(-[rho-0.5].^2)]
%Dibujo del gradiente de T:
quiver(rho,teta,Tx,Ty)
%Fijación de los ejes:
axis([0,4,-2,2])
view(2)

Para comprobar gráficamente que el gradiente de temperatura es ortogonal a las curvas de nivel, creamos el siguiente programa:

Gradiente de la temperatura

x=linspace(0,8,30);
y=linspace(-1,1,30);
rho=sqrt(x.^2+y.^2)
teta=atan(y./[x+0.0000001])
[rho,teta]=meshgrid(x,y)
T=1+[sin(teta).^2].*[exp(-[rho-0.5].^2)]
Tx=-2.*[sin(teta).^2].*[rho-0.5].*[exp(-[rho-0.5].^2)]
Ty=sin(2*teta).*[exp(-[rho-0.5].^2)]
figure(1)
contour(rho,teta,T,100)
figure(2)
quiver(rho,teta,Tx,Ty)
figure(3)
quiver(rho,teta,Tx,Ty)
hold on
contour(rho,teta,T,10,'r')
hold off

10 . Cálculo de la presión media en los puntos del fluido

Para calcular la presión media en los puntos del fluido es necesario aproximar la integral de la presión en todo el fluido y dividirlo por el área total en el canal.

La integral empleada para calcular la presión es una integral de superficie con los límites [0,8],[-1,1], que son los intervalos de los parámetros x e y.

La integral a calcular será la siguiente:

Suponemos [math] p_{1} = 2 [/math] y [math] p_{2} = 1 [/math] y la integral nos queda:

Para resolverla emplearemos el método del Trapecio

, donde y N es el número de puntos entre los que se divide la región a integrar.

Para dicho cálculo, hemos diseñado el siguiente programa en Matlab:

% Número de puntos
N = 200;
a=0; b=8;
h = (b-a)/N;
A = 0;
% Función a integrar
f = @(x) (2*(x-3));
for k = 1:N
x0 = a + (k-1)*h;
x1 = a + k*h;
% Método del trapecio 
A = A + h/2*(f(x0) + f(x1));
end
% Resultado
disp(A)

Con el programa, la integral queda:

Una vez resuelta dicha integral, necesitamos dividir el resultado entre el área del canal, que la obtendremos mediante la siguiente operación:

Sabiendo que se trata de un rectángulo de dimensiones [0,8],[-1,1], el área será:

[math] Área = 8\cdot{2} = 16 [/math]

Una vez obtenidos éstos cálculos, la presión media será:

11 . Caudal del canal

11.1 Caudal del canal

Para calcular el caudal emplearemos la siguiente integral de línea:

siendo t un parámetro y N su vector normal.

Empleando la velocidad y suponiendo [math] p_{1} = 2 [/math], [math] p_{2} = 1 [/math] y [math] \mu = 1 [/math], obtenemos la integral a calcular que nos dará el caudal (sustituyendo 'y' por el parámetro t):

Para calcular dicha intregral volvemos a emplear el método del Trapecio en Matlab:

%Número de puntos
N = 200;
a=-1; b=1;
h=(b-a)/N;
A = 0;
% Función a integrar
f = @(t) ((1-t^2)/2)
for k = 1:N
t0 = a+(k-1)*h;
t1=a+k*h;
% Método del trapecio 
A = A + h/2*(f(t0) + f(t1));
end
% Resultado
disp(A)

Ejecutando este programa, obtenemos que el caudal es [math] 0.6666\frac{M^2}{S} [/math]

11.2 Porcentaje del caudal que pasa por el centro del canal

Por último, calcularemos el porcentaje del caudal que pasa por el centro del canal. Para ello emplearemos la misma integral que para calcular el caudal pero cambiando el intervalo. En éste caso, el intervalo será aquel que tenga el valor medio de y en el centro, siendo -0.5 y 0.5 los extremos.

Modificamos dichos extremos en el mismo programa de Matlab empleado para el caudal

%Número de puntos
N = 200;
a=-0.5; b=0.5;
h=(b-a)/N;
A = 0;
% Función a integrar
f = @(t) ((1-t^2)/2)
for k = 1:N
t0 = a+(k-1)*h;
t1=a+k*h;
% Método del trapecio 
A = A + h/2*(f(t0) + f(t1));
end
% Resultado
disp(A)

Una vez modificado, lo ejecutamos y obtenemos una aproximación de [math] 0.4583\frac{M^2}{S} [/math]. Con este dato, ya podemos obtener el porcentaje pedido. Lo obtendremos mediante la siguiente fórmula:

[math] Caudal(Porcentaje) = \frac{Caudal_{centro}}{Caudal_{total}}\cdot 100 = \frac{0.4583}{0.6666}\cdot 100 ≈ 68.75 [/math]