Visualización de campos escalares y vectoriales en fluidos. (Grupo B6)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Visualización de campos escalares y vectoriales en fluidos. (Grupo B6) |
| Asignatura | Teoría de Campos |
| Curso | 2021-22 |
| Autores | María Magaldi Jurado, Nuria Martínez Ballester, Julia Cabeza Duque |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
En este trabajo vamos a visualizar campos escalares y vectoriales en fluidos. Más específicamente, el fluido a tratar es uno incompresible que se encuentra alrededor de un obstáculo de forma circular, formando así un anillo comprendido entre las circunferencias de radio 1 y 5, con el centro en el origen de coordenadas. Para facilitar el trabajo usaremos coordenadas polares (cilíndricas [math](\rho,\theta)[/math]). Para dibujar que el fluido va por el exterior de la superficie definida consideraremos que [math](\rho,\theta) \in [0,5]*[0,2\pi][/math].
Además, tendremos que la velocidad de las partículas del fluido viene definida por el gradiente de la función potencial [math]\varphi=(\rho+\frac{1}{\rho})sin(\theta)[/math]
Contenido
- 1 º Mallado del fluido
- 2 º Función potencial [math]\varphi[/math] y velocidad de las partículas del fluido
- 3 º Ortogonalidad [math]\vec u[/math] y [math]\vec n[/math]
- 4 º Calculo de [math]\vec u[/math] con [math] \frac{1}{\rho}[/math] despreciable
- 5 º Rotacional y Divergencia
- 6 ºLíneas de corriente
- 7 º Puntos frontera
- 8 º Presión del fluido
- 9 ºEcuación de Navier-Stokes
- 10 º Trayectoria de la línea de corriente
- 11 º Paradoja D'Alembert
- 12 º Curvas de nivel de la presión
1 º Mallado del fluido
El siguiente mallado representa los puntos interiores de la región ocupada por un fluido. Se encuentra en el exterior del círculo unidad. Como ya se ha introducido, nuestro mallado del anillo está comprendido entre los radios 1 y 5 y centro el origen.
El CÓDIGO matlab empleado es:
ro=linspace(1,5,50); %creación de los intervalos
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt); %creación matrices
Mx=Mro.*cos(Mtt); %parametrización de x
My=Mro.*sin(Mtt); %parametrización de y
Mz=zeros(size(Mx)); %parametrización de z (nula al estar en un plano)
subplot(1,1,1);
mesh(Mx,My,Mz) %mallado
axis([-5,5,-5,5]) %superficie definida en [-5,5]*[-5,5]
view(2)
2 º Función potencial [math]\varphi[/math] y velocidad de las partículas del fluido
2.1 º Función potencial
La velocidad de las partículas del fluido viene dada por el gradiente de la función potencial: [math]\varphi=(\rho+\frac{1}{\rho})sin(\theta)[/math] Las curvas de nivel de esta función potencial son las siguientes:
El código matlab empleado es:
ro=linspace(1,5,50);
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt);
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));
f=(Mro+1./Mro).*sin(Mtt); %función potencial
subplot(1,1,1)
contour(Mx,My,f,50); colorbar %curvas de nivel función potencial
2.2 º Campo de velocidades (gradiente de función potencial)
Si dibujamos nuestro campo de velocidades [math]\vec u=\nabla\varphi[/math] en matlab, la representación del mismo con flechas deberá ser ortogonal a las curvas de nivel antes dibujadas. Ya que, de forma geométrica, el gradiente es un vector normal a la curva de nivel de la función de la cual es gradiente en el punto que se esté estudiando.
Hemos pasado nuestra función potencial a cartesianas y trabajado posteriormente con sus derivadas.
Así, con las igualdades: [math]\rho=\sqrt{x^2+y^2}\qquad\theta=arctg(\frac{y}{x})[/math]
Obtenemos la función potencial: [math]\varphi=y+\frac{y}{x^2+y^2}[/math]
El gradiente se define como: [math]\vec u=\nabla\varphi=\frac{∂φ}{∂x}\vec i + \frac{∂φ}{∂j}\vec j [/math]
Realizamos las derivadas parciales:
[math]FX=\frac{\partial\varphi(x,y)}{\partial x}=\frac{-2xy}{(x^2+y^2)^2}\qquad FY=\frac{\partial\varphi(x,y)}{\partial y}=\frac{x^2-y^2}{(x^2+y^2)^2}[/math]
Con estas funciones, podremos dibujar en matlab el campo de velocidades con el siguiente código:
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));
f=(Mro+1./Mro).*sin(Mtt); %función potencial
syms x y %nuevas variables en cartesianas para definir derivadas (campo de velocidades)
fx=inline('-2*(x.*y)./((x.^2+y.^2).^2)','x','y'); %derivadas que definen nuestro campos de velocidades
fy=inline('1+((x.^2-y.^2)./((x.^2+y.^2).^2))','x','y');
d1=fx(Mx,My);
d2=fy(Mx,My);
hold on
subplot(1,1,1)
contour(Mx,My,f,20); colorbar %curvas de nivel función potencial
quiver(Mx,My,d1,d2) %campo de velocidades fluido, ortogonal a curvas de nivel
axis([-5,5,-5,5])
hold offEn nuestro dibujo saldrán las flechas, que representan nuestro campo de velocidades superpuestas con las curvas de nivel de la función potencial. Si ampliamos la imagen podremos observar que estos dos elementos son ortogonales:
3 º Ortogonalidad [math]\vec u[/math] y [math]\vec n[/math]
Como el obstáculo es un círculo situado en el plano xy, sabemos con certeza que entonces el vector normal será el correspondiente al eje z, es decir: [math]\vec n= \vec k[/math] Del apartado anterior conocemos que [math]\vec u=\nabla\varphi[/math].
[math]\vec{u}=sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}[/math]
Conociendo entonces el vector normal y el vector gradiente, podemos calcular su producto escalar. Para comprobar si dos vectores son ortogonales, su producto escalar debe de ser nulo, es por ello que realizamos el producto escalar del gradiente y del vector normal. [math]\nabla\varphi\cdot\vec n[/math] = [math](sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}})\cdot\vec ({e}_{z})[/math] = 0
Como podemos observar, el producto escalar resulta 0, y queda comprobado así que el vector función potencial es perpendicular a su vector normal.
4 º Calculo de [math]\vec u[/math] con [math] \frac{1}{\rho}[/math] despreciable
Al centrarnos en el fluido a una distancia muy lejana, podemos considerar [math] \frac{1}{\rho}[/math] despreciable. En este caso sabremos que [math] u[/math] cambia. La velocidad de las partículas desde los puntos lejanos del fluido viene definida por el gradiente de la función potencial [math]\varphi=\rho\sin(\theta)[/math]
[math]\vec{u}=sin(\theta){\vec{e}_{\rho}}+\frac{1}{\rho}\rho\cos(\theta){\vec{e}_{\theta}}[/math]
En este caso: [math]|\vec{u}| =\sqrt{cos(\theta)^2+sin(\theta)^2}=1[/math]
5 º Rotacional y Divergencia
[math]\nabla\times\vec u=\begin{vmatrix} \vec {e}_{\rho}&\rho\vec {e}_{\theta}&\vec {e}_{z} \\ \frac{\partial}{\partial{\rho}} & \frac{\partial}{\partial{\theta}} & \frac{\partial}{\partial{z}} \\ sin(\theta)[1-\frac{1}{\rho^2}] & \rho\cos(\theta)[\rho +\frac{1}{\rho}] & {0} \end{vmatrix}=[-cos(\theta)(1-\frac{1}{\rho^2})]\vec {e}_{z}+[cos(\theta)(1-\frac{1}{\rho^2})]\vec {e}_{z}=\vec {0}[/math]
[math]\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}({\rho}sin(\theta)(1 -\frac{1}{\rho^2}))+\frac{\partial}{\partial{\theta}}(cos(\theta)(\rho +\frac{1}{\rho})+\frac{\partial}{\partial{z}}(\rho\cdot{0})]=\frac{1}{\rho}[sin(\theta)(1+\frac{1}{\rho^2})-sin(\theta)(1+\frac{1}{\rho^2})]=0 [/math]
En este apartado comprobamos analíticamente que la divergencia y rotacional de [math] u[/math] son nulas. Que la divergencia sea nula indica que en este caso el volumen del fluido no varía, es decir, se mantiene, ni se expande ni se contrae. A esta condición se le denomina incompresibilidad. Además, que el rotacional sea nulo demuestra que las partículas del fluido no giran.
6 ºLíneas de corriente
Las líneas de corriente son aquellas que son tangentes a cada punto del campo [math]\vec u[/math]. Para ello calcularemos el vector [math]\vec v=\vec k\times\vec u[/math]
[math]\vec v=\frac{1}{\rho}\begin{vmatrix} \vec {e}_{\rho}&\rho\vec {e}_{\theta}&\vec {e}_{z} \\ {0} & {0} & {1} \\ sin(\theta)(1-\frac{1}{\rho^2}) & cos(\theta)(1+\frac{1}{\rho^2})\rho & {0} \end{vmatrix}= sin(\theta)(1-\frac{1}{\rho^2})\vec {e}_{\theta} - cos(\theta)(1+\frac{1}{\rho^2})\vec {e}_{\rho} =\vec {grad}\psi[/math]:
[math]\psi=\int - cos(\theta)(1+\frac{1}{\rho^2}) \, d\rho \,\!=-cos(\theta)(\rho-\frac{1}{\rho}) + f(\theta) [/math]:
[math]\frac{\partial\psi}{\partial\theta}=sin(\theta)(\rho-\frac{1}{\rho})+f'(\theta)=sin(\theta)(\rho-\frac{1}{\rho})\qquad f'(\theta)=0\qquad f(\theta)=cte[/math]:
[math]\psi=-cos(\theta)(\rho-\frac{1}{\rho})[/math]:
Del apartado anterior, podemos ver que el rotacional es nulo, y por lo tanto es irrotacional. De esta forma, las líneas de corriente serán tangentes al campo.
Las curvas de nivel de la corriente son las siguientes:
Si superponemos la imagen con el campo de velocidades antes calculado, observamos que efectivamente son tangentes:
El código Matlab empleado es el siguiente:
ro=linspace(1,5,300); %intervalo de ro (1,5)
tt=linspace(0,2*pi,300); %intervalo de teta (0,2pi)
[Mro,Mtt]=meshgrid(ro,tt);
Mx=Mro.*cos(Mtt);
My=Mro.*sin(Mtt);
Mz=zeros(size(Mx));
syms x y
fi=inline('-x+x./(x.^2+y.^2)','x','y'); %función que define la corriente de u
f2=fi(Mx,My);
fx=inline('-2*(x.*y)./((x.^2+y.^2).^2)','x','y'); %derivadas de funcion potencial que definen campo de velocidades del fluido
fy=inline('1+((x.^2-y.^2)./((x.^2+y.^2).^2))','x','y');
d1=fx(Mx,My);
d2=fy(Mx,My);
hold on
subplot(1,1,1)
contour(Mx,My,f2,100) %curvas de nivel de la corriente
quiver(Mx,My,d1,d2) %campo de velocidades fluido, tangente a lineas de corriente
hold off
7 º Puntos frontera
Al tener definida la superficie del anillo con los radios 1 y 5, buscamos donde están situados los puntos de frontera. Con la función potencial inicial [math]\varphi=(\rho+\frac{1}{\rho})sin(\theta)[/math] y con [math] u[/math] = [math](sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}) [/math].
Calculamos el modulo del vector en función de [math]\theta[/math] para
[math]\rho=1\qquad [/math] [math] u[/math] = [math] 2cos(\theta){\vec{e}_{\theta}} [/math]
[math]|\vec u|=\sqrt{4cos(\theta)^2} = 2cos(\theta) = 2|cos(\theta)|[/math]
Como podemos observar, el módulo de la velocidad es máximo cuando cos[math]\theta[/math] alcanza sus valores máximos (en 0 y en [math]\pi[/math], donde alcanza en valor absoluto el valor de 1).
Así, los puntos de velocidad máxima son [math]\theta=0[/math] y [math]\theta={\pi}[/math].
Los puntos de velocidad mínima o puntos de remanso son [math]\theta=\frac{\pi}{2}[/math] y [math]\theta=\frac{3\pi}{2}[/math].
8 º Presión del fluido
En dinámica de fluidos, el principio de Bernoulli, también denominado ecuación de Bernoulli, describe el comportamiento de un fluido moviéndose a lo largo de una línea de corriente. Relacionando las siguientes variables:
[math]p =[/math]presión estática a la que está sometido el fluido debida a las moléculas que lo rodean.
[math]d =[/math]densidad del fluido.
[math]\vec { u } =[/math]velocidad de flujo del fluido. :
[math]\frac { 1 }{ 2 } d { \left| \vec { u } \right| }^{ 2 }+ { p } =cte[/math]
Para que esta expresión pueda aplicarse, el fluido debe ser incompresible, no viscoso y fluir en régimen estacionario (la velocidad en un punto determinado no varía con el tiempo). Tomaremos nuestro fluido como tal.
Suponiendo que nuestro fluido efectivamente cumple la ecuación de Bernouilli, que posee una densidad constante de d=2 y tomando la cte de la ecuación como cte=10, calcularemos la presión de nuestro fluido.
Sustituyendo los valores mencionados tenemos que:
[math] p = 10- |\vec{u}|^2 [/math]
Calculamos el módulo al cuadrado de nuestro campo de velocidades:
[math] |\vec{u}|^2 = (1-\frac{1}{ρ^2})^2 sen^2(θ) + (1+\frac{1}{ρ^2})^2 cos^2(θ) [/math]
Así la función que define la presión del fluido es:
[math] p = 10- (1-\frac{1}{ρ^2})^2 sen^2(θ) + (1+\frac{1}{ρ^2})^2 cos^2(θ) [/math]
Para la representación de la presión empleamos el siguiente código:
ro=linspace(1,5,50);
tt=linspace(0,2*pi,50);
[Mro,Mtt]=meshgrid(ro,tt);
x=Mro.*cos(Mtt);
y=Mro.*sin(Mtt);
P=10-(((1-1./Mro.^2).^2).*(sin(Mtt).^2)+((1+1./Mro.^2).^2).*(cos(Mtt).^2)); %función presión
surf(x,y,P); colorbar;
axis([-5,5,-5,5])
view(2)
Pmax=max(max(P)) %presión máx y mín
Pmin=min(min(P))Matlab nos devuelve la siguiente imagen y los valores máximos y mínimos de la presión (Presión máx: 9.9959 Presión mínima: 6)
Recordemos ahora la imagen del campo de velocidades (ahora sin ampliar) que habíamos obtenido anteriormente con el código matlab de el apartado 2.2:
Si comparamos nuestra imagen de la presión del fluido (Fig.PRESIÓN DEL FLUIDO) con la del campo de velocidades (Fig.CAMPO DE VELOCIDADES) podemos observar que los puntos que cuentan con módulo de velocidad máxima ([math]\theta=0[/math] y [math]\theta={\pi}[/math]) son aquellos con presiones mínimas. Esto es un resultado coherente que verifica la ecuación de Bernoulli, a mayor velocidad, menores presiones.
La ecuación de Navier-Stokes describe el movimiento de los fluidos a partir de su campo de velocidades y de presiones:
[math] d(\vec{u}·\nabla)\vec{u} + \nabla p = \mu \Delta \vec{u} [/math]
Los fluidos incompresibles tienen la propiedad de ocupar siempre el mismo volumen, de esta forma la viscosidad será nula: [math] \mu =0 [/math] y la densidad será nula también: [math]d = 2[/math]. Si sustituimos en la ecuación general, llegamos a la conclusión de que para un fluido incompresible:
[math] 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 [/math]
Si consideramos los campos [math] (\vec{u},p) [/math] hallados anteriormente:
[math] (\vec{u}·\nabla)\vec{u} = u_{1} \frac{\partial u_{1}}{\partial \rho} \vec{e}_{\rho} + u_{1} \frac{\partial u_{1}}{\partial \theta} \vec{e}_{\theta}+ u_{2} \frac{\partial u_{2}}{\partial \rho} \vec{e}_{\rho}+ u_{2} \frac{\partial u_{2}}{\partial \theta} \vec{e}_{\theta} = [/math] :
[math] (\vec{u}·\nabla)\vec{u}=(sen(\theta)(1-\frac{1}{\rho^2})[\frac{2}{\rho^3}sin(\theta)\vec{e}_{\rho}+(1-\frac{1}{\rho^2})cos(\theta)\vec{e}_{\theta}])+((1+\frac{1}{\rho^2})cos(\theta)[-\frac{2}{\rho^3}cos(\theta)\vec{e}_{\rho}-(1+\frac{1}{\rho^2})sin(\theta)\vec{e}_{\theta}])=[/math]
[math] =(\vec{u}·\nabla)\vec{u}=[\frac{2}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)-\frac{2}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[(1-\frac{1}{\rho^2})^2sin(\theta)cos(\theta)-(1+\frac{1}{\rho^2})^2sin(\theta)cos(\theta)]\vec{e}_{\theta} [/math]
Lo aplicamos a nuestro caso con la densidad determinada:
[math] 2(\vec{u}·\nabla)\vec{u}=[\frac{4}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)- \frac{4}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[(1-\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)-(1+\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)]\vec{e}_{\theta} [/math]
De la misma forma:
[math] \nabla p = \frac{\partial p}{\partial \rho} \vec{e}_{\rho} + \frac{\partial p}{\partial \theta} \vec{e}_{\theta} = [/math] :
[math] \nabla p = -[(2(1-\frac{1}{\rho^2})\frac{2}{\rho^3}sin^2(\theta))+(2(1+\frac{1}{\rho^2})(-\frac{2}{\rho^3})cos(\theta))]\vec{e}_{\rho} - [((1-\frac{1}{\rho^2})^2)2sin(\theta)cos(\theta)+((1+\frac{1}{\rho^2})^2)2(-sin(\theta))cos(\theta)]\vec{e}_{\theta}=[/math]
[math] = [-\frac{4}{\rho^3}(1-\frac{1}{\rho^2})sin^2(\theta)+ \frac{4}{\rho^3}(1+\frac{1}{\rho^2})cos^2(\theta)]\vec{e}_{\rho}+[-(1-\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)+(1+\frac{1}{\rho^2}^2)2sin(\theta)cos(\theta)]\vec{e}_{\theta} [/math]
Habiendo sustituido todos los valores en la ecuación de Navier-Stokes, coparando los resultados vemos directamente que se cumple [math] 2(\vec{u}·\nabla)\vec{u} + \nabla p = 0 [/math]
quedando así demostrado que se cumple la ecuación.
10 º Trayectoria de la línea de corriente
Si fuéramos una partícula del fluido con el que estamos trabajando, seguiríamos la trayectoria de una línea de corriente, es decir empezaríamos con su trayectoria original y a medida que se va acercando al obstáculo recorreríamos el contorno de este, y cuando lo hubiésemos pasado, volveríamos a la dirección inicial. En nuestro caso, como el obstáculo tiene un contorno circular el fluido empieza recto, luego se curva para evitar el obstáculo y finalmente vuelve a ir en la dirección inicial. En relación a la variación de la velocidad y presión podemos comprobar con las gráficas, que en los puntos donde la velocidad es mínima las presiones son máximas. En los puntos que ya hemos calculado anteriormente donde la velocidad es 0, [math]\theta=\frac{\pi}{2}[/math] y [math]\theta=\frac{3\pi}{2}[/math], sabemos que la presión es máxima. Con esto podemos deducir que la presión y velocidad de estas partículas son inversamente proporcionales.
11 º Paradoja D'Alembert
[math]\vec{u}=sin(\theta)(1-\frac{1}{\rho^2}){\vec{e}_{\rho}}+cos(\theta)(1+\frac{1}{\rho^2}){\vec{e}_{\theta}}[/math]
Para [math]\rho=1[/math]:
[math]\vec u=2cos(\theta)\vec {e}_{\theta}[/math]:
Según el Teorema de Kutta-Joukowski la fuerza que ejerce el fluido sobre el obstáculo, en este caso el círculo de radio 1, es proporcional a la circulación. Por lo tanto:
[math]\int_{0}^{2\pi} \vec u r'(\theta)\,d\theta =\int_{0}^{2\pi} 2cos(\theta)\vec{e}_{\theta}\vec{e}_{\theta}\,d\theta =\int_{0}^{2\pi} 2cos(\theta)\,d\theta=2[sin(0)-sin(2\pi)]=0 [/math]:
El fluido no ejerce ninguna fuerza sobre el obstáculo, dando lugar a la paradoja D'Alembert, que concluye que la fuerza resultante sobre el cuerpo es cero, contradiciéndose con la observación. Esto demuestra también que la circulación del campo vectorial a lo largo de la circunferencia se anula. Para poder ilustrar esta fuerza se necesitan tomar soluciones de modelos de fluidos con una viscosidad positiva.
12 º Curvas de nivel de la presión
Observando las curvas de nivel de la presión podemos entender mejor las conclusiones a las que hemos ido llegando a lo largo del trabajo. En esta gráfica podemos observar los puntos de presión máxima [math]\theta=\frac{\pi}{2}[/math] y [math]\theta=\frac{3\pi}{2}[/math] y los puntos de presión mínima [math]\theta=0[/math] y [math]\theta=\pi[/math]. Además, cabe destacar que, en las zonas en las que las curvas de nivel están mas juntas, hay mayor variación de la presión.
Código Matlab:
contour(x,y,P,100) % Curvas de nivel de la presión
axis([-5,5,-5,5]) % Región a dibujar
