Visualización de un campo escalar y vectorial en un fluido (Grupo G7)

Trabajo realizado por estudiantes
Título	Visualización de un campo escalar y vectorial en un fluido. Grupo G7
Asignatura	Teoría de Campos
Curso	2014-15
Autores	Álvaro Baeza Cabrero, Daniel Fojo Berlana, Pablo Carrasco del Olmo, Jorge Juan Fernández Díaz
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

En este artículo se va a realizar el estudio del flujo de un fluido incompresible a través de un canal.

1.1 Concepto de fluido incompresible

Un fluido incompresible es un tipo de medio continuo cuyas moléculas tienen poca cohesión entre sí y cuya densidad permanece aproximadamente constante a lo largo de todo el fluido, condición necesaria para considerar un fluido incompresible. Los fluidos están conformados por los líquidos y los gases, siendo los segundos mucho menos viscosos (casi fluidos ideales). Los líquidos toman la forma del recipiente que los aloja, manteniendo su propio volumen (ausencia de memoria de forma).Esto se debe a su separación molecular en la que los fluidos no poseen una forma definida por tanto no se puede calcular su volumen o densidad a simple vista, para esto se introduce el fluido en un recipiente en el cual toma su forma y así podemos calcular su volumen y densidad, esto facilita su estudio.

Fluido

Fluido en recipiente

Algunas de las propiedades que más se suelen tener en cuenta en un fluido serían:

1. Presión
2. Temperatura
3. Densidad
4. Compresibilidad
5. Viscosidad

Ahora pasamos al estudio del fluido de este canal, en algunas de estas propiedades y de otras propiedades.

2 Campo de velocidades y presiones

Supongamos que el canal por el que se mueve el fluido tiene paredes rectas de dimensiones [0,4]x[0,1]. Vamos a representar el canal con el fluido dentro de él, mediante Matlab, cuyo código sería:

Representación del canal y del fluido

% Se crea un vector llamado x con valores entre 0 y 4 y con una equidistancia de 0.1

x=0:0.1:4;

% Creamos un vector y con valores entre 0 y 1 y de equidistancia 0.1

y=0:0.1:1; 

% Realizamos el mallado de los vectores x e y

[xx,yy]=meshgrid(x,y); 

% Se Dibuja el mallado

mesh(xx,yy,xx*0) 

% Los ejes los fijamos en [0,4]x[-1,2]

axis([0,4,-1,2]) 

% Representamos el campo en el plano OXY

view(2)

Este fluido en particular, tiene los siguientes campos:

[math] p(x,y)=p_{1}+(p_{2}-p_{1})(x-1) [/math]

Sería un campo escalar, es el campo de presiones del fluido que estamos estudiando.

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

Es un campo vectorial, representa el campo de velocidades de este fluido.

Siendo [math] p_{1}[/math] la presión en los puntos cuya abcisa es [math]x=1[/math], [math] p_{2}[/math] la presión en los puntos [math] x=2[/math], y [math] \mu[/math] el coeficiente de viscosidad del fluido. Podemos ver que la velocidad en las paredes del canal es nula, en [math]y=0[/math] e [math]y=1[/math] ya que: [math] \vec u (x,0)= \frac{0(1-0)(p_{1}-p_{2})}{(2\mu)}\vec i = \vec 0 [/math]: [math] \vec u (x,1)= \frac{1(1-1)(p_{1}-p_{2})}{(2\mu)}\vec i = \vec 0 [/math]

2.1 Ecuación de Navier-Stokes

A continuación se comprueba que el fluido sigue un régimen estacionario, para ello se tiene que verificar la ecuación de Navier-Stokes:

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

Empezamos calculando [math] \nabla \vec u [/math] que se haría del siguiente modo:

[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{(1-2y)(p_{1}-p_{2})}{(2\mu)}\vec j [/math]

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

Esto se debe a que el gradiente del campo es ortogonal al propio campo, por lo que su producto escalar sería 0. Ahora vamos a calcular [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]

Ahora pasamos al segundo miembro de la ecuación, calculando [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 = \frac{-(p_{1}-p_{2})}{(\mu)} \vec i [/math]

Como vemos se verifica la ecuación de Navier-Stokes, por lo que el fluido sigue un régimen estacionario, como queríamos demostrar.

2.2 Verificación de la condición de incompresibilidad del fluido

Para que un fluido sea incompresible se tiene que verificar que [math] \nabla \cdot \vec u = 0[/math], es decir la divergencia del campo de velocidades tiene que ser nula, pasamos a comprobarlo:

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

La condición se cumple y en consecuencia de ello, el fluido es incompresible.

2.3 Representación de campos de velocidades y presiones

Para representar el campo de velocidades y presiones vamos a suponer que [math] p_{1}=2[/math], [math] p_{2}=1[/math] y [math] \mu=1[/math],siendo el campo de presiones y velocidades:

[math] p(x,y)=3-x [/math]

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

Para representar los campos hemos utilizado matlab cuyos códigos son respectivamente:

Campo de presiones

% Vector x con valores entre 0 y 4, y equidistancia 0.1

x=0:0.1:4;  

% Vector y con valores entre 0 y 1, y equidistancia 0.1

y=0:0.1:1; 

% Mallado de la superficie

[xx,yy]=meshgrid(x,y); 
figure (1)

% Campo de presiones particularizado

p=3-xx  

% Representación del campo escalar de presiones
      
surf(xx,yy,p)      

% Fijación de los ejes en [0,4]x[-1,2]

axis([0,4,-1,2])
 
% Visualización del campo en el plano OXY
 
view(2)

Campo de velocidades

% Vector x con valores entre 0 y 1, y equidistancia 0.1 

x=0:0.1:4;   

% Vector y con valores entre 0 y 1, y equidistancia 0.1       

y=0:0.1:1;   

% Matrices de las coordenadas x e y

[xx,yy]=meshgrid(x,y); 
figure (2)

% Campo de velocidades en x

ux=(1/2).*(yy-(yy.^2)); 

% Campo de velocidades en y

uy=0.*xx;
               
% Dibujamos el campo vectorial de velocidades

quiver(xx,yy,ux,uy)

% Fijación de los ejes en [0,4]x[-1,2]

axis([0,4,-1,2])

% Visualización del campo en el plano OXY  

view(2)

En el campo de presiones podemos ver que en las zonas más cálidas (colores como rojo, naranja, amarillo y verde), la presión es más alta y en las zonas más frías (colores como verde,celeste y azul), la presión es más baja, esto significa que la presión baja a medida que el fluido va avanzando en el canal.

En el campo de velocidades, apreciamos que todos los vectores apuntan a la misma dirección, paralelo al eje x. También que la velocidad es mayor en el centro del canal que en las paredes porque los vectores tienen mayor módulo ahí, esto ocurre porque no hay barreras en el centro que pueda provocar una disminución en la velocidad.

2.4 Velocidad máxima del fluido

La velocidad será máxima en el punto cuya derivada de su módulo sea cero y cuya derivada segunda en el punto de máxima velocidad tenga un valor inferior a cero. Es decir:

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

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

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

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

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

En general será máximo o mínimo en función de los valores de [math] p_{1} [/math], [math] p_{2} [/math] y [math] \mu [/math]. En el caso de que [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 tanto, en ese caso, la velocidad sería máxima en todos los puntos cuya ordenada sea [math] y = \frac{1}{2} [/math]. Este método se puede utilizar debido a que el campo depende de una sola variable, si no fuera así habría que calcular la matriz hessiana, que es el gradiente de un gradiente: [math] H = \nabla (\nabla \vec u) = \begin{bmatrix} \frac{\partial^2 \vec u}{\partial x^2} & \frac{\partial^2 \vec u}{\partial x \partial y} \\\frac{\partial^2 \vec u}{\partial x \partial y} & \frac{\partial^2 \vec u}{\partial y^2} \end{bmatrix} [/math]

3 Líneas de corriente

Las líneas de corriente siguen una dirección ortogonal al campo de velocidades, se hallaría del siguiente modo:

[math] \vec{v}=\vec{k}\times \vec{u} =\frac{y(1-y)(p_1-p_2)}{2\mu}\vec{j} [/math]

Ese campo que hemos hallado por ser ortogonal al anterior eso significa que [math] \vec{v} = \nabla \psi [/math] en el que [math] \psi [/math] es la función de corriente. Vamos a proceder a calcular esta función:

[math] \vec v = \frac{ \partial \psi}{\partial x} \vec i + \frac{ \partial \psi}{\partial y} \vec j [/math]

[math] \frac{ \partial \psi}{\partial x} = v_{1} = 0 [/math]

[math] \frac{ \partial \psi}{\partial y} = v_{2} = \frac{y(1-y)(p_1-p_2)}{2\mu} \implies \psi = \int\frac{y(1-y)(p_1-p_2)}{2\mu}dy = \frac{(p_1-p_2)}{2\mu}(\frac{y^2}{2}-\frac{y^3}{3})+ f(x) [/math]

[math] \frac{ \partial \psi}{\partial x} = 0 + f’(x) = 0 \implies f’(x) = 0 \implies f(x) = K [/math]

Siendo K una constante. Entonces las líneas de corriente siguen la siguiente función:

[math] \psi = \frac{(p_1-p_2)}{2\mu}(\frac{y^2}{2}-\frac{y^3}{3}) + K [/math]

Para dibujar las líneas de corriente vamos a suponer que [math] p_{1} = 2 [/math], [math] p_{2} = 1 [/math], [math] \mu = 1 [/math] y [math] K = 0 [/math], siendo su función potencial [math]\psi =\frac{1}{2}(\frac{y^2}{2}-\frac{y^3}{3})[/math]. Para representarlo hemos utilizado matlab con el siguiente código:

Líneas de corriente del campo de velocidades

% Vector x con valores comprendidos entre 0 y 4 y de equidistancia 0.1

x=0:0.1:4;

% Vector y con valores comprendidos entre 0 y 1 y de equidistancia 0.1
           
y=0:0.1:1;             

% Matrices de las coordenadas x e y para el mallado

[xx,yy]=meshgrid(x,y); 

% función del potencial escalar

lineas_de_corriente= (1/2).* ((yy.^2/2)-(yy.^3/3) ); 

% Se dibuja las líneas de corriente en el plano OXY

contour(xx,yy,lineas_de_corrientes) 

% Fijación de los ejes en [0,4]x[-1,2]

axis([0,4,-1,2])          

% Visualización de las líneas de corriente en el plano OXY

view(2)

Ahora vamos a demostrar que las líneas de corriente del campo de velocidades [math] \vec u [/math] es lo mismo que calcular la función potencial del campo [math] \vec v [/math]:

Llamemos ε a la función de corriente de [math] \vec u [/math]

[math] \vec{v} = \vec{k}\times \vec{u}= v_1\vec{i} + v_2\vec{j} = -u_2\vec{i} +u_1\vec{j} [/math]

Siendo [math] u_1 = \frac{ \partial ε}{\partial y} [/math], [math] u_2 = -\frac{ \partial ε}{\partial x} [/math],[math] v_1 = \frac{ \partial \psi}{\partial x} [/math] y [math] v_2 = \frac{ \partial \psi}{\partial y} [/math]. Como [math] v_1 = - u_2 [/math] y [math] v_2 = u_1 [/math], entonces:

[math] \frac{ \partial \psi}{\partial x} = -(-\frac{ \partial ε}{\partial x}) [/math]

[math] \frac{ \partial \psi}{\partial y} = \frac{ \partial ε}{\partial y}\implies \psi = ε + K [/math]

Esto demuestra que es lo mismo las líneas de corriente de [math] \vec u [/math] que la función potencial de [math] \vec v [/math].

3.1 Campo irrotacional

Como se demostró antes [math] \nabla \cdot \vec u = \frac{\partial u_1}{\partial x} + \frac{\partial u_2}{\partial y} = 0[/math], vamos a ver si esto implica que el campo [math] \vec v [/math] sea irrotacional: [math]\nabla\times\vec{v}= det\begin{pmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ v_1 & v_2 & v_3 \end{pmatrix}= det\begin{pmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ -u_2 & u_1 & 0 \end{pmatrix}=-\frac{\partial}{\partial z}(u_1)\vec{i} + \frac{\partial}{\partial z}(-u_2)\vec{j} + (\frac{\partial}{\partial x}(u_1) -\frac{\partial}{\partial y}(-u_2))\vec{k} = 0\vec {i} + 0\vec{j} + (\frac{\partial}{\partial x}(u_1) -\frac{\partial}{\partial y}(-u_2))\vec{k} = 0\vec {i} + 0\vec{j} + (\nabla \cdot \vec u)\vec{k} = \vec 0 [/math] Como queríamos demostrar el campo [math] \vec v [/math] es irrotacional por ser su rotación igual a cero.

4 Rotacional del campo de velocidades

Vamos a calcular el rotacional del campo [math] \vec u [/math], después se calcula su módulo y debido a que depende de una sola variable el campo se puede hallar el punto de máxima rotacional derivando el módulo del rotacional e igualándola a cero: [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-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] [math]|\nabla\times\vec{u}| = \frac{(1-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]

Esta igualdad implica que no existen máximos ni mínimos relativos. Vamos a ver en las paredes si son máximos o mínimos absolutos, por ser extremos de un intervalo cerrado y acotado:

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

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

En el caso particular [math] p_{1} = 2 [/math], [math] p_{2} = 1 [/math] y [math] \mu = 1 [/math]:

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

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

Por lo tanto, los puntos de mayor rotacional se encuentran en las paredes del canal, esto se debe a que la velocidad es mayor en la zona central del canal que en las paredes, por lo que se ve razonable estos resultados. A continuación vamos a proceder a representar el módulo del rotacional del campo [math] \vec u [/math], que sería la función [math]|\nabla\times\vec{u}| = \frac{1-2y}{2} [/math]. Utilizamos matlab para ello con el siguiente código:

Representación gráfica del módulo del rotacional del campo de velocidades

% Vector x con valores entre 0 y 4 y equidistancia 0.1

x=0:0.1:4;            

% Vector y con valores entre 0 y 1 y equidistancia 0.1

y=0:0.1:1;           

% Matrices de las coordenadas x e y 

[xx,yy]=meshgrid(x,y); 

% Módulo del rotacional del campo de velocidades

modulo_rotacional=abs (1/2.*(1-2.*yy)) 

% Dibujo de la superficie   

surf(xx,yy,modulo_rotacional)     

% Fijación de los ejes en [0,4]x[-1,2]

axis([0,4,-1,2]) 

% Visualización del módulo del rotacional en el plano OXY

view (2)

5 Campo de temperaturas

5.1 Representación del campo de temperaturas

Ahora vamos a hablar del campo de temperaturas, un campo escalar que te indica la temperatura en cada punto del fluido. En este canal el campo de temperaturas T(x,y) sigue la siguiente función: [math] T(x,y) = e^{-(x-1)^{2}+y^{2}} [/math] Vamos a pasar a representar el campo de temperaturas mediante matlab con el siguiente código:

Campo de temperaturas

% Vector x con valores entre 0 y 4 de equidistancia 0.1

x= 0:0.1:4; 

 % Vector y con valores entre 0 y 1 y equidistancia 0.1
 
y= 0:0.1:1; 

% Matrices de las coordenadas x e y para el mallado

[xx,yy]=meshgrid(x,y); 

% Campo de temperaturas

T =exp(-(xx-1).^2 +(yy.^2)) 

% Dibujo del campo escalar

surf(xx,yy,T)         

% Fijación de los ejes en [0,4]x[-1,2]
axis([0,4,-1,2])
% Visualización del campo de temperaturas en el plano OXY
view(2)

Como podemos ver en el gráfico que la temperatura máxima se encuentra alrededor de [math] x = 1 [/math] en la pared superior y si vamos avanzando hacia la izquierda del canal, la temperatura del fluido va disminuyendo.

También vamos a representar las curvas de nivel del campo de temperaturas, que se utiliza el siguiente código de matlab:

Curvas de nivel del campo de temperaturas

% Vector x con valores entre 0 y 4 con una equidistancia de 0.1

x= 0:0.1:4;      

% Vector y con valores entre 0 y 1 y de equidistancia 0.1

y= 0:0.1:1;      

% Matrices de las coordenadas x e y para el mallado

[xx,yy]=meshgrid(x,y); 

% Campo de temperaturas

T=exp(-(xx-1).^2 +(yy.^2)); 

% Dibujo de las curvas de nivel 

contour (xx,yy,T)    

% Fijación de los ejes en [0,4]x[-1,2]
axis([0,4,-1,2])
% Visualización de las curvas de nivel y del gradiente de temperaturas en el plano OXY
view(2)

De igual modo podemos ver por el color de las curvas que alrededor de [math] x = 1 [/math] y en la pared superior la temperatura es más alta por tener un color más cálido ahí (rojo).

5.2 Representación del gradiente de temperaturas

Ahora vamos a hallar el gradiente del campo de temperaturas. Por ser un campo escalar, se calcularía del siguiente modo: [math] \nabla {T(x,y)} = \frac{\partial{T}}{\partial x}\vec{i} + \frac{\partial{T}}{\partial y}\vec{j} = -2(x-1)\cdot{e^{-(x-1)^{2}+y^{2}}}\vec{i} + 2y\cdot{e^{-(x-1)^{2}+y^{2}}}\vec{j} [/math]

Vamos a representarlo en un gráfico con el siguiente código de matlab:

Gradiente del campo de temperaturas

% Vector x con valores entre 0 y 4 de equidistancia 0.1

x=0:0.1:4;      

% Vector y con valores entre 0 y 1 y equidistancia 0.1

y=0:0.1:1;      

% Matrices de las coordenadas x e y para el mallado

[xx,yy]=meshgrid(x,y); 

% Primer componente del gradiente

Ti=-2.*(xx-1).*exp(-(xx-1).^2 +(yy.^2)); 

% Segunda componente del gradiente
Tj= 2.*yy.*exp(-(xx-1).^2 +(yy.^2)); 
% Dibujo del gradiente del campo de temperaturas
quiver(xx,yy,fi,fj) 
% Fijación de los ejes en [0,4]x[-1,2]
axis([0,4,-1,2])
% Visualización del gradiente en el plano OXY
view(2)

A continuación vamos a representar las curvas de nivel del campo de temperaturas y el gradiente del campo en un mismo gráfico y comprobaremos si son ortogonales. Para ello utilizamos matlab introduciendo el siguiente código:

Curvas de nivel y el gradiente del campo de temperaturas

% Vector x con valores entre 0 y 4 con una equidistancia de 0.1

x= 0:0.1:4;       

% Vector y con valores entre 0 y 1 y de equidistancia 0.1

y= 0:0.1:1;       

% Matrices de las coordenadas x e y para el mallado

[xx,yy]=meshgrid(x,y);
 
% Campo de temperaturas

T=exp(-(xx-1).^2 +(yy.^2)); 

% Dibujo de las curvas de nivel

contour (xx,yy,T)     

% Para representar las curvas de nivel y el gradiente de temperatura al mismo tiempo

hold on 

% Primer componente del gradiente

Ti= -2.*(xx-1).*exp(-(xx-1).^2 +(yy.^2));  

% Segunda componente del campo gradiente

Tj= 2.*yy.*exp(-(xx-1).^2 +(yy.^2)); 

% Dibujo del gradiente del campo de temperaturas   
quiver (xx,yy,Ti,Tj)
% Fijación de los ejes en [0,4]x[-1,2]
axis([0,4,-1,2]) 
% Visualización de las curvas de nivel y del gradiente de temperaturas en el plano OXY
view(2)

Podemos ver que efectivamente el gradiente es ortogonal al campo de temperaturas, ya que cada uno de los vectores que se componen el gradiente siguen una dirección perpendicular a la que siguen las curvas de nivel del campo de temperaturas.

6 Integrales numéricas

6.1 Presión media

Vamos a calcular la presión media de los puntos del fluido en el canal. Para ello vamos a calcular la integral de la presión en todo el fluido y dividirlo por el área del canal en el que se encuentra. Empecemos por la integral de la presión que será una integral de superficie siendo los parámetros x e y y sus intervalos respectivamente son [0,4] y [0,1], por lo que la integral que habrá que calcular será la siguiente: [math] \int_S p(x,y) \; dS = \int_0^4(\int_0^1 p_{1}+(p_{2}-p_{1})(x-1) \; dy) \;dx = \int_0^4 p_{1}+(p_{2}-p_{1})(x-1) \; dx [/math]

Supongamos que [math] p_{1} = 2 [/math], [math] p_{2} = 1 [/math], entonces:

[math] \int_0^4 (3-x) \; dx [/math]

Esta integral la vamos a aproximar por la regla del trapecio en el que:

[math] \int_0^4 p(x) \; dx \sim h\frac12 p(x_0)+h\sum_{i=1}^{N-1}p(x_i)+h\frac12 p(x_N) [/math]

Siendo [math] h = \frac{4-0}{N} [/math] y siendo N el número de trapecios que vas a dividir la región a integrar, cuantos más sean, la aproximación va a ser mejor. Para hacer este cálculo hemos utilizado matlab con el siguiente código:

N = 200;% Número de puntos
a=0; b=4;% Extremos del intervalo
h = (b-a)/N;% La longitud de cada trapecio
A = 0;
f = @(x) (3-x);% Función a integrar
for k = 1:N
x0 = a + (k-1)*h;% Valor inicial de la función en una cantidad de trapecios entre 1 y 200
x1 = a + k*h;% Valor final de la función en una cantidad de trapecios entre 1 y 200
A = A + h/2*(f(x0) + f(x1));% Regla del trapecio 
end
disp(A)% Resultado

Utilizando este método, la integral nos da el siguiente resultado:

[math] \int_0^4 (3-x) \; dx = 4 [/math]

Ahora hay que dividirlo entre el área del canal que se calcularía del siguiente modo:

[math] \int_S \; dS = Área [/math]

Pero como el área es un rectángulo de dimensiones [0,4]x[0,1], el área sería igual a:

[math] Área = 4\cdot{1} = 4 [/math]

Por lo que la presión media será igual a:

[math] P_{media} = \frac{\int_S p(x,y) \; dS}{\int_S \; dS} = \frac{4}{4} = 1 [/math]

Por lo que la presión media nos daría un resultado de 1.

6.2 Caudal del canal

El caudal se calcularía con una integral de línea del siguiente modo:

[math] Caudal = \int_γ \vec {u(t)}\cdot \vec N \; ds [/math]

Siendo t un parámetro y [math] \vec N [/math] su vector normal. En este caso, [math] y = t [/math] con un intervalo entre 0 y 1. El vector normal se calcula a partir del vector tangente del siguiente modo:

[math] \vec t = \frac{d\vec r}{dt} = A\vec i + B\vec j [/math]

Siendo [math] \vec r [/math] el vector posición:

[math] \vec N = -B\vec i + A\vec j [/math]

El vector posición sería con un valor de x cualquiera por ejemplo supongamos que fuera [math] x = 2 [/math], entonces el vector posición sería: [math] \vec r = 2\vec i + t\vec j \implies \vec N = -\vec i [/math]

Así que la integral a calcular será la siguiente en el caso [math] p_{1} = 2 [/math], [math] p_{2} = 1 [/math] y [math] \mu = 1 [/math]:

[math] \int_0^1 (\frac{-t(1-t)}{2}) \; dt [/math]

Dará un valor negativo, por lo que hemos dado el sentido de la normal opuesta a la que tiene, para cambiar el sentido sólo cambiamos el signo. Para el cálculo de esta integral vamos a utilizar el método del trapecio, con el siguiente código matlab:

N = 200;% Número de puntos
a=0; b=1;% Extremos del intervalo
h=(b-a)/N;% La longitud de cada trapecio
A = 0;
f = @(t) ((t*(1-t))/2)% Función a integrar
for k = 1:N
t0 = a + (k-1)*h;% Valor inicial de la función en una cantidad de trapecios entre 1 y 200
t1 = a + k*h;% Valor final de la función en una cantidad de trapecios entre 1 y 200
A = A + h/2*(f(t0) + f(t1));% Regla del trapecio 
end
disp(A)% Resultado

Con este código el caudal nos daría una aproximación de [math]0.083331\frac{M^2}{S}[/math], las unidades en la que está expresado el caudal se debe a que se aplica sobre una sección.

6.2.1 Porcentaje de caudal

Para concluir con el estudio de este fluido se va a calcular el porcentaje que pasa por el centro del canal. El caudal que pasa por el centro se calculará con la misma integral que se aplicó en el caudal total, pero su intervalo será entre 0.25 y 0.75, teniendo el valor medio de y, en el centro del intevalo. El código de matlab utilizado sería el siguiente:

N=200;% Número de puntos
a=0.25; b=0.75;% Extremos del intervalo
h=(b-a)/N;% La longitud de cada trapecio
A = 0;
f = @(t) ((t*(1-t))/2);% Función a integrar
for k = 1:N
t0 = a + (k-1)*h;% Valor inicial de la función en una cantidad de trapecios entre 1 y 200
t1 = a + k*h;% Valor final de la función en una cantidad de trapecios entre 1 y 200
A = A + h/2*(f(t0) + f(t1));% Regla del trapecio 
end
disp(A)% Resultado

El caudal nos daría una aproximación de [math] 0.057291\frac{M^2}{S} [/math] Ahora vamos a calcular el porcentaje de caudal que pasa por el centro: [math] Caudal(Porcentaje) = \frac{Caudal_{centro}}{Caudal_{total}}\cdot 100 = \frac{0.057291}{0.083331}\cdot 100 ≈ 68.75 [/math] Viendo el resultado que nos da vemos que más de la mitad del caudal total pasa por el centro, creemos que se debe a que la velocidad es mayor en el centro que en las paredes.