Trabajo De Campos Grupo 7

Trabajo realizado por estudiantes
Título	Visualización de un campo escalar y vectorial en un fluido. Grupo 7-C
Asignatura	Teoría de Campos
Curso	2019-20
Autores	Natalia Peña Salas, Stella Bárbara Peña Martínez, Alejandro Vegazo Luengo
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 . Mallado que representa los puntos ocupados por un fluido

Dibujaremos 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]. El código para crear el programa aparece a continuación, junto con su representación gráfica:

%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)

2 . Ecuación de Navier-Stokes

Verificamos 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.

3 . 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 . 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].

6 . 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])

7 . 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)

8 . 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)

9 . 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á:

10 . Caudal del canal

10.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]

10.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]