function tcm12
%TCM12 resuelve el problema 12 de transporte de cantidad de movimiento.
%
%   El objetivo del problema será obtener el perfil de velocidad y de
%   densidad de flujo de cantidad de movimiento para el caso de flujo
%   de un fludio entre dos placas paralelas tal y como se a descrito en el 
%   problema 12. Los datos que se conocerán son los siguientes:
%
%      - Perdida de presión por metro de longitud: -dP/dx=6870 (N/m2)/m
%      - Distancia entre placas: 0.02 m
%      - Viscosidad dinámica del fluido (aceite): mu=0.18 kg/(m·s)
%      - Densidad del fluido (aceite): ro=850 kg/m3
%      - Aceleración de la gravedad: g=9.81 m/s2
%      - Ángulo de inclinación respecto a la horizontal: alfa=30º=pi/6
%
%Aplicaremos la herramienta de resolución (solver) BVP4C para resolver
%problemas de ecuaciones diferenciales ordinarias (Ordinary Differential
%Equations, ODEs) con condiciones de valores límite (Boundary Value 
%Problem, BVP). 

% NOTA: A nivel interno de cálculo, y dado que la estructura "solinit" que
% Matlab utiliza en BVP4C maneja como variables "x,y", tomaremos 
% y=velocidad y x=z.

dPX = -6870;
a = -0.01;
b = 0.01;
mu = 0.18;
ro = 850;
g = 9.81;
alfa = pi/6;
% Definimos los datos del problema

function dydx = odefun_tcm12(x,y)
%ODEFUN_TCM12  Define el sistema de ODEs del problema TCM12. 
%   Matlab sólo trabaja con ecuaciones diferenciales de primer orden, por
%   lo que si aparecen ecuaciones diferenciales de orden superior es
%   preciso transformarlas en un sistema de ecuaciones difereciales de
%   primer orden. Para ello se realizará un cambio de variables llamando
%   y(1)=y, y(2)=y', y(3)=y'', y(4)=y''', y así sucesivamente; en nuestro 
%   caso como tenemos una ecuación diferencial de orden dos es preciso 
%   definir sólo y(1)=y, y(2)=y'. El sistema de ecuaciones diferenciales de
%   primer orden del problema TCM12 quedaría definido con la siguiente 
%   expresión:
dydx =  [ y(2)
         (dPX-ro*g*sin(alfa))/mu ];
end
     
 function res = bcfun_tcm12(ya,yb)
%BCFUN_TCM12  Define las condiciones límite del problema TCM12.
%   RES = BCFUN_TCM12 devuelve el vector columna RES con los residuos
%   resultantes de las soluciones aproximadas YA e YB en los extremos del
%   intervalo definido por las condiciones límite [a b]. El problema BVP 
%   está resuelto cuando RES=0 (es decir cuando se anula la diferencia 
%   entre la solución aproximada y el valor de la solución en los puntos 
%   correspondientes a las condiciones límite).
res = [ ya(1)
        yb(1)];
end
    
function yinit = yinit_tcm12(x)
%YINIT_TCM12 Solución aproximada (guess) para el problema TCM12.
%   Para resolver problemas BVP es preciso introducir una aproximación a la
%   solución mediante un vector columna YINIT. La aproximación será bien
%   una aproximación al valor de la función y al valor de la primera
%   derivada de la función en la condición límite xa, o bien la expresión 
%   de una función que aproxime el comportamiento de la solución buscada; 
%   en el problema TCM12 hemos optado por esta última opción sabiendo que 
%   el perfil de velocidades es parabólico y presenta su máximo en el punto 
%   medio entre las condiciones límite, por lo que una ecuación polinómica 
%   de segundo orden con su máximo en x=0 m, es decir "-x^2" será adecuada.
%   Tomando este criterio el primer elemento del vector columna YINIT será 
%   esta función y el segundo elemento su primera derivada:
yinit = [ -x^2
          -2*x ];
end

solinit_tcm12 = bvpinit(linspace(a,b,10),@yinit_tcm12);
%BVPINIT Es una función auxiliar para crear la estructura inicial de la
%solución aproximada.
%   Una vez tenemos YINIT_TCM12 Matlab precisa crear una matriz inicial de
%   valores aproximados de la solución para diferentes valores de x
%   dentro del intervalo definido por las condiciones límite. Esa matriz se
%   construye utilizando la función auxiliar BVPINIT, cuya solución
%   almacenaremos en SOLINIT_TCM12. La sintaxis de BVPINIT es la siguiente:
%       
%       bvpinit(x,yinit)
%
%           - x: valores de x para la matriz inicial. Aquí es donde se
%                definen que valores de x comprendidos entre las 
%                condiciones límite [a b] ensayaremos para calcular la 
%                solución inicial aproximada. En nuestro caso nos hemos 
%                basado en el comando LINSPACE para crear 10 puntos 
%                linealmente espaciados entre A y B.
%           -yinit: aproximación del valor de la solución en los valores de 
%                x de la matriz inicial. En nuestro caso la hemos definido 
%                en la función previa YINIT_TCM12. 

sol = bvp4c(@odefun_tcm12,@bcfun_tcm12,solinit_tcm12);
%BVP4C es la herramienta de resolucion de Matlab para BVPs.
%   La sintaxis de BVP4C es la siguiente:
%
%       bvp4c(odefun,bcfun,solinit)
%
%           - odefun: sistema de ecuaciones diferenciales de primer orden
%                a solucionar; en nuestro caso ODEFUN_TCM12.
%           - bcfun: evaluación de los residuos en las condiciones límite;
%                en nuestro caso BCFUN_TCM12.
%           - solinit: estructa de la solución inicial aproximada; en
%                nuestro caso SOLINIT_TCM12.
%
%   La solución generada se almacena en la estructura SOL, donde:
%           - sol.x: son los valores de x seleccionados por BVP4C para
%                evaluar la solución.
%           - sol.y: la solución del problema en los puntos sol.x.
%           - sol.yp: evaluación de la primera derivada de la solución, y',
%                en los puntos sol.x.
%           - sol.solver: nombre de la herramienta utilizada para resolver
%                el problema, en nuestro caso 'bvp4c'.
%
xint = linspace(a,b,101);
[Sxint, SPxint] = deval(sol,xint);
TAU = -mu*SPxint(1,:);
%PERFIL DE VELOCIDADES: Para construir el perfil de velocidades 
%   utilizaremos la función auxiliar DEVAL que nos devuelve el valor de la
%   solución del problema para los valores de x que deseemos. En nuestro 
%   caso dividiremos el intervalo entre A y B en subintervalos de tamaño
%   0.0002 m mediante el comando LINSPACE, guardando el valor de cada 
%   punto en XINT. Almacenaremos los valores de la solución en la variable
%   SXINT (los valores se alamcenan en la primera fila de SXINT).
%PERFIL DE DENSIDAD DE FLUJO DE CANTIDAD DE MOVIMIENTO, TAU: La orden DEVAL
%   también nos permite obtener el valor de la primera derivada del
%   polinomio que interpolaría la solución entre los puntos de la
%   solución, resumiendo la primera derivada de la solución, la cual se
%   almacena en la variable SPXINT(los valores se alamcenan en la primera 
%   fila de SPXINT). Si multiplicamos SPXINT por -MU obtendremos el perfil 
%   de TAU para nuestro caso particular. 
plot(Sxint(1,:),xint);
grid
title('Perfil Velocidades')
ylabel('z (m)')
xlabel('Vx (m/s)')
figure(2)
plot(TAU,xint,'r')
grid
title('Perfil Densidad de Flujo de Cantidad de Movimiento')
ylabel('z (m)')
xlabel('TAUzx (N/m2)')
global RESULTADOS12;
RESULTADOS12=[ xint' Sxint(1,:)' TAU'];
RESULTADOS12_VAR=[ '      z(m)       ' 'Vx(m/s)   ' 'TAUzx (N/m2)'];
format short e
disp(RESULTADOS12_VAR)
disp(RESULTADOS12)
%REPRESENTACIÓN Y OBTENCIÓN DE LOS RESULTADOS: La secuencia de comandos
%anterior permite la representación de gráficos de los perfiles, la
%obtención de los resultados en una tabla y el almacenamiento de estos en
%la variable global RESULTADOS12.
end
%
%EJECUCIÓN DE TCM12: Ejecutar en la ventana de comandos:
%
%                global RESULTADOS12; tcm12
%
% Obtenemos los resultados en una tabla en la ventana de comandos, los 
% gráficos de los perfiles y los resultados almacenados en la variable
% global RESULTADOS12 que ahora tenemos disponible en el espacio de trabajo 
% por si queremos hacer algún cálculo posterior con ellos.
%
% 000000000000000000000000000000000
% ---------------------------------
% Juan F. Martínez Gallegos
% Dpto. Ingeniería Química
% UGR
% ---------------------------------
% 000000000000000000000000000000000
