Diferencia entre revisiones de «Difusión de una sustancia contaminante (Grupo 6)»
(→Tiempo en alcanzar la solución estacionaria) |
(→Tiempo en alcanzar la solución estacionaria) |
||
| Línea 719: | Línea 719: | ||
===Tiempo en alcanzar la solución estacionaria=== | ===Tiempo en alcanzar la solución estacionaria=== | ||
El estado estacionario es el valor de concentración u = 3 como podemos comprobar en el gráfico y obteniendo los valores en cualquier punto de x para un valor de tiempo alto. Para encontrar el tiempo en que se alcanza este valor consideramos un error del 5%, que corresponderá a 0,15. Por tanto buscaremos el tiempo a partir del cual todos los puntos de la barra tienen su concentración comprendida en el intervalo [2,85-3,15]. Ayudándonos por la matriz solución y sabiendo que los extremos de la barra corresponden a los valores máximo y mínimo de la barra en todo instante (ver gráfico), tendiendo ambos asintóticamente al 3 desde los valores 2 y 4, basta con comprobar cuando los valores en los extremos llegan a entrar en el intervalo. Dado que 3 es el valor medio entre los valores iniciales en ambos extremos 2 y 4, la aproximación se produce de forma simétrica y en el punto de columna 364 de la matriz solución encontramos los valores dentro del intervalo pedido. Pidiendo a matlab el tiempo en que se produce este elemento con el comando t(364) obtenemos el valor de tiempo 13.77233 segundos, que es el tiempo que tarda en alcanzar el estado estacionario con error del 5%. | El estado estacionario es el valor de concentración u = 3 como podemos comprobar en el gráfico y obteniendo los valores en cualquier punto de x para un valor de tiempo alto. Para encontrar el tiempo en que se alcanza este valor consideramos un error del 5%, que corresponderá a 0,15. Por tanto buscaremos el tiempo a partir del cual todos los puntos de la barra tienen su concentración comprendida en el intervalo [2,85-3,15]. Ayudándonos por la matriz solución y sabiendo que los extremos de la barra corresponden a los valores máximo y mínimo de la barra en todo instante (ver gráfico), tendiendo ambos asintóticamente al 3 desde los valores 2 y 4, basta con comprobar cuando los valores en los extremos llegan a entrar en el intervalo. Dado que 3 es el valor medio entre los valores iniciales en ambos extremos 2 y 4, la aproximación se produce de forma simétrica y en el punto de columna 364 de la matriz solución encontramos los valores dentro del intervalo pedido. Pidiendo a matlab el tiempo en que se produce este elemento con el comando t(364) obtenemos el valor de tiempo 13.77233 segundos, que es el tiempo que tarda en alcanzar el estado estacionario con error del 5%. | ||
| + | |||
| + | |||
La solución cambia ligeramente si ∆x se divide entre diez ya que la aproximación al resultado real será mejor. Pero la variación es muy pequeña dado que ya teníamos un ∆x suficientemente pequeño para la resolución. Lo encontramos en la columna 3695 y nos da un valor de 13.8523 segundos. | La solución cambia ligeramente si ∆x se divide entre diez ya que la aproximación al resultado real será mejor. Pero la variación es muy pequeña dado que ya teníamos un ∆x suficientemente pequeño para la resolución. Lo encontramos en la columna 3695 y nos da un valor de 13.8523 segundos. | ||
Revisión del 20:13 28 abr 2017
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Difusión de una sustancia contaminante. Grupo B-6 |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2016-17 |
| Autores | Daniel Pacheco Sánchez 917 Oscar Lázaro González 993 Alonso Herranz Hudson 1043 Manuel Bécares Martín 1077 Pablo Morales Santón 1177 Dariusz Adam Pabian 1187 |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
El estudio del presente trabajo es el análisis matemático de fenómeno de difusión de una sustancia contaminante a lo largo de un tubo y conceptos relacionados con el mismo, como el empleo de limpiadores.
Contenido
1 Interpretación del fenómeno físico y planteamiento del problema matemático
El fenómeno físico de la difusión de sustancias se rige por las leyes de Fick, que reciben su nombre de Adolf Fick, quien las derivó en 1855. Estas leyes determinan que cuando se da una situación en la que existe una variación de concentración de una sustancia se produce un flujo de partículas que tiende a provocar que la disolución se homogeneice hasta una situación con concentración uniforme a lo largo de todo el medio. La primera ley de Fick establece que el flujo de difusión del contaminante [math]F[/math] es proporcional a la variación de concentración [math]u[/math], suponiendo un estado estacionario:
[math]F=-D\cdot\frac{\partial u}{\partial x}[/math]
donde [math]D[/math] es el coeficiente de difusión medido en [math]\frac{m^2}{s}[/math].
El caso a estudiar es el de un tubo largo en el cual se encuentra una solución compuesta por dos sustancias de las cuales una de ella es un contaminante. Como ya hemos mencionado para enunciar la ley de Fick, denominaremos [math]u[/math] a la concentración de contaminante en cada posición del tubo, medida en [math]\frac{mol}{m^3}[/math]. La longitud del tubo [math]L[/math] será de 8 metros y ocupará el intervalo [math]x \in (0,L)[/math] (se orienta el tubo según el eje x) y su sección transversal será constante a lo largo de toda su longitud. Se supone que la concentración es la misma en cualquier punto de la sección transversal del tubo. Por tanto, la concentración dependerá únicamente de dos variables: [math]u = u(x,t)[/math]. En los extremos se ha colocado un aislante que impide que se produzca flujo hacia el exterior del tubo (por tanto este será nulo).
Se define el flujo de contaminante [math]F(x,t)[/math] como la cantidad del mismo que atraviesa una sección transversal por unidad de tiempo y área, medido en número de moles. Dado que la tubería no permite el flujo a través de su pared y sus extremos están también aislados, el flujo solo se producirá en la dirección [math]x[/math] y dentro de la longitud de la tubería.
Para la obtención de la ecuación que define el proceso de difusión del contaminante tomaremos como leyes físicas que rigen el proceso: el principio de conservación de la masa y la primera ley de Fick, ya enunciada.
El principio de conservación de la masa nos permite determinar que la variación de la cantidad de contaminante por unidad de tiempo en un volumen infinitesimal de tubo es igual a la suma del flujo de contaminante a través de los extremos del volumen por unidad de tiempo, mas la concentración de contaminante que se genera o se pierde en el interior del volumen por unidad de tiempo. Esta última pérdida o ganancia se supondrá igual a cero en el caso estudiado.
Para expresarlo matemáticamente, se calcula primeramente cual es la cantidad de sustancia contaminante en un volumen infinitesimal:
Cantidad de contaminante: [math] A \cdot u(x,t) \cdot Δx [/math] , donde: [math]A[/math] es la sección del tubo medida en [math]m^2[/math] y [math]u(x,t)[/math] esla concentración del contaminante en [math]\frac{mol}{m^3}[/math] .
Si derivamos respecto respecto al tiempo: [math] A \cdot u_t(x,t) \cdot Δx [/math] , que representa la variación de la cantidad de contaminante respecto del tiempo multiplicado por el volumen.
Suponemos [math] Δx \gt 0 [/math] y que la concentración del contaminante en un instante [math] t [/math] es menor en [math] x+Δx [/math] que en [math] x [/math], lo que implica que [math] u(x+Δx,t) - u(x,t) \lt 0 [/math], y al ser [math] Δx [/math] muy pequeño, se tiene que [math]u_x(x,t)\lt0[/math] y el flujo de difusión es positivo según el eje [math]x[/math].
El flujo de sustancia contaminante en el volumen considerado será: [math] F(x,t) \cdot A - F(x+Δx,t) \cdot A ± f(x,t)[/math]
Dado que no se produce una pérdida o ganancia en la concentración del contaminante: [math]f(x,t) = 0[/math].
Por tanto, según el principio de conservación de la masa: [math] A \cdot u_t(x,t) \cdot Δx = F(x,t) \cdot A - F(x+Δx,t) \cdot A[/math]
Dividiendo por [math]Δx[/math] y haciendo que tienda a 0:
[math]u_t(x,t)= \frac {F(x,t)-F(x+Δx, t)}{Δx}[/math]
Aplicando la primera ley de Fick:
[math]u_t(x,t)=-\frac{\partial }{\partial x}(-D \cdot u_x(x,t)) = D \cdot u_{xx}(x,t) [/math]
Se obtiene la ecuación diferencial que rige la difusión de la sustancia contaminante a lo largo del tubo:
[math]u_t(x,t)- D \cdot u_{xx}(x,t)= 0 [/math]
2 Problema propuesto
La ecuación diferencial tiene infinitas soluciones. Para obtener una única solución, se necesitan dos condiciones de frontera y una condición inicial. En el problema que se nos proporciona, tenemos las siguientes condiciones:
- Dado que los extremos están aislados del exterior, no permiten que se produzca flujo a través de ellos, por tanto: [math]u_x(0,t) = u_x(8,t) = 0[/math] .
- En el instante inicial se verifican las siguientes concentraciones: [math] u(x,0)=u_0=\left\{\begin{matrix}2, x≤4\\4, x\gt 4\end{matrix}\right. [/math]
Entonces, el problema a resolver será el siguiente:
- [math] (P)\left\{\begin{matrix}u_t-u_{xx}=0 ,x\in\ (0,8), t\gt0\\u_x(0,t)=0, t\gt0\\u_x(8,t)=0, t\gt0\\u(x,0)=u_0, x\in\ (0,8)\end{matrix}\right. [/math]
2.1 Resolución numérica del problema
Para resolver el problema planteado, se empleará el método numérico de diferencias finitas, el cual busca expresar el problema continuo de partida como un sistema de ecuaciones diferenciales de primer orden.
Para ello, se consideran las siguientes aproximaciones de las derivadas espaciales en x:
- [math]u_{x}(x,t)\simeq\frac{u(x_{n+1},t)-u(x_{n-1},t)}{2h}[/math]
- [math]u_{xx}(x,t)\simeq\frac{u(x_{n-1},t)-2u(x_n,t)+u(x_{n+1},t)}{h^2}[/math]
Aplicando la ecuación diferencial de nuestro sistema en los nodos interiores de la longitud del tubo, obtenemos un sistema de [math]N-1[/math] ([math]N[/math] número de subintervalos de la variable [math]x[/math]) ecuaciones del tipo:
- [math]u_{t}(x,t) + \frac{-u(x_{n-1},t)+2u(x_n,t)-u(x_{n+1},t)}{h^2}=0[/math]
Para las condiciones de contorno, se aplica la aproximación de la primera derivada de [math]u[/math]:
- [math]u'_{0}(t) + \frac{-u_{-1}(t)+2u_{0}(t)-u_{1}(t)}{h^2}=0[/math]
- [math]u_{x}(0,t) = 0 = \frac{u_{1}(t)-u_{-1}(t)}{2h}[/math]
Despejando [math]u_{-1}[/math] en la segunda ecuación y sustituyendo en la primera, se obtiene la ecuación relativa al extremo izquierdo del tubo:
- [math]u'_{0}(t) + \frac{2u_{0}(t)-2u_{1}(t)}{h^2}=0[/math]
Aplicando el mismo procedimiento a la condición en el otro extremo del tubo, se obtiene su ecuación asociada:
- [math]u'_{N+1}(t) + \frac{2u_{N+1}(t)-2u_{N}(t)}{h^2}=0[/math]
Con estas ecuaciones, tenemos el siguiente sistema de [math]N+1[/math] ecuaciones:
- [math] \begin{cases}u'_{0}(t) + \frac{2u_{0}(t)-2u_{1}(t)}{h^2}=0\\u'_{n}(x,t) + \frac{-u(x_{n-1},t)+2u(x_n,t)-u(x_{n+1},t)}{h^2}=0\\u'_{N+1}(t) + \frac{2u_{N+1}(t)-2u_{N}(t)}{h^2}=0\end{cases} [/math]
El cual se puede simplificar como:
- [math]\begin{cases} U'(t)+KU(t)=F(t)=0\\ U(0)=U^0 \end{cases} [/math]
Y este puede resolverse por métodos como: Euler explícito, Euler implícito, trapecio o Heun.
2.1.1 Resolución por el método del trapecio
Para la resolución por los diferentes métodos, se emplearán las siguientes condiciones: [math]\triangle t= \triangle x/4[/math] en [math]t\in\ [0,5][/math], con [math]\triangle x=0.1[/math], aplicadas al sistema ya planteado:
- [math] (P)\left\{\begin{matrix}u_t-u_{xx}=0 ,x\in\ (0,8), t\gt0\\u_x(0,t)=0, t\gt0\\u_x(8,t)=0, t\gt0\\u(x,0)=u_0, x\in\ (0,8)\end{matrix}\right. [/math]
El código MATLAB empleado para su resolución es el siguiente:
clear, close all
%% Datos del problema
a = 0; %Extremo derecho del tubo
b = 8; %Extremo izquierdo del tubo
t0 = 0;
tM = 7; %Tiempo estudiado
D = 1;
%% Funciones
% g : función del término independiente (x vector, t escalar)
g = @(x,t) x*0;
% ca : funcion f(t) de la condicion de contorno en x=a (t vector)
ca = @(t) t*0;
% cb : funcion f(t) de la condicion de contorno en x=b (t vector)
cb = @(t) t*0;
% u0 : forma analitica f(x) funcion valor inicial (x vector)
u0 = @(x) (2.*(x<=4))+(4.*(x>4));
%% Discretización del vector de espacio
incrX = 0.15;
N = round((b-a)/incrX);
x = linspace(a,b,N+1);
x = x';
%% Planteamiento del sistema a resolver U' = -K*U + G
% Matriz K
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2) = -2;
K(N+1,N) = -2;
K = K*(D/incrX^2);
% Condición inicial
U0 = u0(x);
%% Resolucion del sistema como PVI
% Discretización del vector de tiempos
incrT = incrX/4;
M = round((tM-t0)/incrT);
t = linspace(t0,tM,M+1);
%Inicializacion de la matriz de soluciones
sol = zeros(N+1,M+1);
sol(:,1)= U0;
% Método de trapecio
for i=1:M
Gi = g(x,t(i)); % Gi es el término independiente de la ecuación en el elemento n de la iteración (nulo en este caso)
Gi(1) = Gi(1)-2*D*ca(t(i))/incrX;
Gi(end) = Gi(end)+2*D*cb(t(i))/incrX;
Gi2 = g(x,t(i+1)); % Gi2 es el término independiente de la ecuación en el elemento n+1 de la iteración (nulo en este caso)
Gi2(1) = Gi2(1)-2*D*ca(t(i+1))/incrX;
Gi2(end) = Gi2(end)+2*D*cb(t(i+1))/incrX;
sol(:,i+1) = (eye(size(K))+incrT/2*K)\(sol(:,i)+incrT/2*(Gi2-K*sol(:,i)+Gi));
end
%% Representacion grafica
[Mt,Mx]= meshgrid(t,x);
mesh(Mt,Mx,sol)
title('Distribución de las concentraciones de contaminante en el tubo en función del tiempo');
xlabel('Tiempo (s)'); ylabel('Coordenada x (m)'); zlabel('Concentración (mol/(m^2*s))');
Con el cual se obtiene la siguiente gráfica en la que se representan la concentración de contaminante en función del espacio y del tiempo:
En la gráfica correspondiente se puede observar un doble análisis, en función de la proximidad al punto medio del tubo, podemos observar como el punto x=4, supone el limite entre ambos lados, lo que supone que los puntos mas cercamos a él, van a sufrir una variación con un mayor gradiente, lo cuál tiene su lógica física, ya que son los primeros puntos en sufrir la variación y tienen que aumentar lo suficientemente rápido para transmitir ésta variación a los puntos siguientes, hasta los llamados puntos frontera. El otro análisis que se puede hacer, es en función del tiempo, en el que se observa que la gráfica se acerca al principio de forma muy rápida, y se va estabilizando poco a poco hasta que de forma casi tangencial alcanza el valor limite que coincide con el valor medio, que corresponde a su vez a la posición de equilibrio de concentración de contaminante entre ambas mitades del tubo, es decir, que todas las secciones del tubo tengan el mismo valor de concentración. La variación en el gradiente será mayor, contra mayor es la diferencia de concentración de contaminante entre los puntos [math]x=0[/math], y [math]x=8[/math], al igual que contra más se aproximan ambas mitades al valor medio, es cuando varía de forma más lenta. El valor se alcanzará cuándo la masa total de contaminante contenido en el tubo, se distribuya de manera constante a lo largo del mismo.
2.1.2 Resolución por el método de Euler explícito
El código MATLAB empleado para su resolución es el siguiente:
clear, close all
%% Datos del problema
a = 0; %Extremo derecho del tubo
b = 8; %Extremo izquierdo del tubo
t0 = 0;
tM = 7; %Tiempo estudiado
D = 1;
%% Funciones
% g : función del término independiente (x vector, t escalar)
g = @(x,t) x*0;
% ca : funcion f(t) de la condicion de contorno en x=a (t vector)
ca = @(t) t*0;
% cb : funcion f(t) de la condicion de contorno en x=b (t vector)
cb = @(t) t*0;
% u0 : forma analitica f(x) funcion valor inicial (x vector)
u0 = @(x) (2.*(x<=4))+(4.*(x>4));
%% Discretización del vector de espacio
incrX = 0.15;
N = round((b-a)/incrX);
x = linspace(a,b,N+1);
x = x';
%% Planteamiento del sistema a resolver U' = -K*U + G
% Matriz K
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2) = -2;
K(N+1,N) = -2;
K = K*(D/incrX^2);
% Condición inicial
U0 = u0(x);
%% Resolucion del sistema como PVI
% Discretización del vector de tiempos
incrT = incrX/4;
M = round((tM-t0)/incrT);
t = linspace(t0,tM,M+1);
%Inicializacion de la matriz de soluciones
sol = zeros(N+1,M+1);
sol(:,1)= U0;
% Método de Euler explicito
for i=1:M
Gi = g(x,t(i)); % Gi es el término independiente de la ecuación en el elemento n de la iteración (nulo en este caso)
Gi(1) = Gi(1)-2*D*ca(t(i))/incrX;
Gi(end) = Gi(end)+2*D*cb(t(i))/incrX;
sol(:,i+1) = sol(:,i)+incrT*(-K*sol(:,i)+Gi);
end
%% Representacion grafica
[Mt,Mx]= meshgrid(t,x);
mesh(Mt,Mx,sol)
title('Distribución de las concentraciones de contaminante en el tubo en función del tiempo');
xlabel('Tiempo (s)'); ylabel('Coordenada x (m)'); zlabel('Concentración (mol/(m^2*s))');
Con el cual se obtiene la siguiente gráfica en la que se representan la concentración de contaminante en función del espacio y del tiempo:
Como podemos observar en el gráfico, el resultado obtenido por el método de Euler muestra su inestabilidad. Esta situación se debe a que el tamaño de paso temporal ∆t = ∆x/4 es mayor que el valor límite que determina si un método explícito es estable o no (∆t = 0.5 * ∆x^2).
2.1.3 Resolución por el método de Euler implícito
El código MATLAB empleado para su resolución es el siguiente:
clear, close all
%% Datos del problema
a = 0; %Extremo derecho del tubo
b = 8; %Extremo izquierdo del tubo
t0 = 0;
tM = 7; %Tiempo estudiado
D = 1;
%% Funciones
% g : función del término independiente (x vector, t escalar)
g = @(x,t) x*0;
% ca : funcion f(t) de la condicion de contorno en x=a (t vector)
ca = @(t) t*0;
% cb : funcion f(t) de la condicion de contorno en x=b (t vector)
cb = @(t) t*0;
% u0 : forma analitica f(x) funcion valor inicial (x vector)
u0 = @(x) (2.*(x<=4))+(4.*(x>4));
%% Discretización del vector de espacio
incrX = 0.15;
N = round((b-a)/incrX);
x = linspace(a,b,N+1);
x = x';
%% Planteamiento del sistema a resolver U' = -K*U + G
% Matriz K
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2) = -2;
K(N+1,N) = -2;
K = K*(D/incrX^2);
% Condición inicial
U0 = u0(x);
%% Resolucion del sistema como PVI
% Discretización del vector de tiempos
incrT = incrX/4;
M = round((tM-t0)/incrT);
t = linspace(t0,tM,M+1);
%Inicializacion de la matriz de soluciones
sol = zeros(N+1,M+1);
sol(:,1)= U0;
% Método de Euler implícito
for i=1:M
Gi = g(x,t(i)); % Gi es el término independiente de la ecuación en el elemento n de la iteración (nulo en este caso)
Gi(1) = Gi(1)-2*D*ca(t(i))/incrX;
Gi(end) = Gi(end)+2*D*cb(t(i))/incrX;
Gi2 = g(x,t(i+1)); % Gi2 es el término independiente de la ecuación en el elemento n+1 de la iteración (nulo en este caso)
Gi2(1) = Gi2(1)-2*D*ca(t(i+1))/incrX;
Gi2(end) = Gi2(end)+2*D*cb(t(i+1))/incrX;
sol(:,i+1) = (eye(size(K))+incrT*K)\(sol(:,i)+incrT*Gi2);
end
%% Representacion grafica
[Mt,Mx]= meshgrid(t,x);
mesh(Mt,Mx,sol)
title('Distribución de las concentraciones de contaminante en el tubo en función del tiempo');
xlabel('Tiempo (s)'); ylabel('Coordenada x (m)'); zlabel('Concentración (mol/(m^2*s))');
Con el cual se obtiene la siguiente gráfica en la que se representan la concentración de contaminante en función del espacio y del tiempo:
En la nueva resolución por el método de Euler implícito se observa con respecto a la resolución por el método del trapecio que es más exacta en los primeros instantes del tiempo. En estos existía una serie de picos en la concentración en puntos cercanos al medio del tubo que ahora desaparecen.
2.1.4 Resolución por el método de Heun
El código MATLAB empleado para su resolución es el siguiente:
clear, close all
%% Datos del problema
a = 0; %Extremo derecho del tubo
b = 8; %Extremo izquierdo del tubo
t0 = 0;
tM = 7; %Tiempo estudiado
D = 1;
%% Funciones
% g : función del término independiente (x vector, t escalar)
g = @(x,t) x*0;
% ca : funcion f(t) de la condicion de contorno en x=a (t vector)
ca = @(t) t*0;
% cb : funcion f(t) de la condicion de contorno en x=b (t vector)
cb = @(t) t*0;
% u0 : forma analitica f(x) funcion valor inicial (x vector)
u0 = @(x) (2.*(x<=4))+(4.*(x>4));
%% Discretización del vector de espacio
incrX = 0.15;
N = round((b-a)/incrX);
x = linspace(a,b,N+1);
x = x';
%% Planteamiento del sistema a resolver U' = -K*U + G
% Matriz K
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2) = -2;
K(N+1,N) = -2;
K = K*(D/incrX^2);
% Condición inicial
U0 = u0(x);
%% Resolucion del sistema como PVI
% Discretización del vector de tiempos
incrT = incrX/4;
M = round((tM-t0)/incrT);
t = linspace(t0,tM,M+1);
%Inicializacion de la matriz de soluciones
sol = zeros(N+1,M+1);
sol(:,1)= U0;
% Método de Heun
for i=1:M
Gi1 = g(x,t(i)); % Gi1 es el vector G en t(i)
Gi1(1) = Gi1(1)-2*D*ca(t(i))/incrX;
Gi1(end) = Gi1(end)+2*D*cb(t(i))/incrX;
Gi2 = g(x,t(i)+incrT); %Gi2 es el vector G en t(i)+k
Gi2(1) = Gi2(1)-2*D*ca(t(i)+incrT)/incrX;
Gi2(end) = Gi2(end)+2*D*cb(t(i)+incrT)/incrX;
K1 = -K*sol(:,i)+Gi1;
K2 = -K*(sol(:,i)+K1*incrT)+Gi2;
sol(:,i+1) = sol(:,i)+incrT/2*(K1+K2);
end
%% Representacion grafica
[Mt,Mx]= meshgrid(t,x);
figure(1)
mesh(Mt,Mx,sol)
title('Distribución de las concentraciones de contaminante en el tubo en función del tiempo');
xlabel('Tiempo (s)'); ylabel('Coordenada x (m)'); zlabel('Concentración (mol/(m^2*s))');
Con el cual se obtiene la siguiente gráfica en la que se representan la concentración de contaminante en función del espacio y del tiempo:
Como en el caso del método de Euler explícito, este método es también inestable, pero con unos picos de mayor magnitud ya que el gráfico muesta el eje de la concentración en 10^226
2.2 Comparativa de métodos empleados
A continuación se analizan los errores como medio para comparar los diferentes métodos. Los gráficos de errores obtenidos para realizar la comparativa se han obtenido con el siguiente código MATLAB:
%% Datos del problema
a = 0; %Extremo derecho del tubo
b = 8; %Extremo izquierdo del tubo
t0 = 0;
tM = 7; %Tiempo estudiado
D = 1;
%% Funciones
% g : función del término independiente (x vector, t escalar)
g = @(x,t) x*0;
% ca : funcion f(t) de la condicion de contorno en x=a (t vector)
ca = @(t) t*0;
% cb : funcion f(t) de la condicion de contorno en x=b (t vector)
cb = @(t) t*0;
% u0 : forma analitica f(x) funcion valor inicial (x vector)
u0 = @(x) (2.*(x<=4))+(4.*(x>4));
%% Discretización del vector de espacio
incrX = 0.15;
N = round((b-a)/incrX);
x = linspace(a,b,N+1);
x = x';
%% Planteamiento del sistema a resolver U' = -K*U + G
% Matriz K
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2) = -2;
K(N+1,N) = -2;
K = K*(D/incrX^2);
% Condición inicial
U0 = u0(x);
%% Resolucion del sistema como PVI
% Discretización del vector de tiempos
incrT = incrX/4;
M = round((tM-t0)/incrT);
t = linspace(t0,tM,M+1);
%Inicializacion de la matriz de soluciones
soltrap = zeros(N+1,M+1); soleuler = zeros(N+1,M+1); soleulerimp = zeros(N+1,M+1); solheun = zeros(N+1,M+1);
soltrap(:,1)= U0; soleuler(:,1)= U0; soleulerimp(:,1)= U0; solheun(:,1)= U0;
% Método de trapecio
for i=1:M
Gi = g(x,t(i)); % Gi es el término independiente de la ecuación en el elemento n de la iteración (nulo en este caso)
Gi(1) = Gi(1)-2*D*ca(t(i))/incrX;
Gi(end) = Gi(end)+2*D*cb(t(i))/incrX;
Gi2 = g(x,t(i+1)); % Gi2 es el término independiente de la ecuación en el elemento n+1 de la iteración (nulo en este caso)
Gi2(1) = Gi2(1)-2*D*ca(t(i+1))/incrX;
Gi2(end) = Gi2(end)+2*D*cb(t(i+1))/incrX;
soltrap(:,i+1) = (eye(size(K))+incrT/2*K)\(soltrap(:,i)+incrT/2*(Gi2-K*soltrap(:,i)+Gi));
end
% Método de Euler
clear Gi Gi2
for i=1:M
Gi = g(x,t(i)); % Gi es el término independiente de la ecuación en el elemento n de la iteración (nulo en este caso)
Gi(1) = Gi(1)-2*D*ca(t(i))/incrX;
Gi(end) = Gi(end)+2*D*cb(t(i))/incrX;
soleuler(:,i+1) = soleuler(:,i)+incrT*(-K*soleuler(:,i)+Gi);
end
% Método de Euler implícito
clear Gi Gi2
for i=1:M
Gi = g(x,t(i)); % Gi es el término independiente de la ecuación en el elemento n de la iteración (nulo en este caso)
Gi(1) = Gi(1)-2*D*ca(t(i))/incrX;
Gi(end) = Gi(end)+2*D*cb(t(i))/incrX;
Gi2 = g(x,t(i+1)); % Gi2 es el término independiente de la ecuación en el elemento n+1 de la iteración (nulo en este caso)
Gi2(1) = Gi2(1)-2*D*ca(t(i+1))/incrX;
Gi2(end) = Gi2(end)+2*D*cb(t(i+1))/incrX;
soleulerimp(:,i+1) = (eye(size(K))+incrT*K)\(soleulerimp(:,i)+incrT*Gi2);
end
% Método de Heun
clear Gi Gi2
for i=1:M
Gi1 = g(x,t(i)); % Gi1 es el vector G en t(i)
Gi1(1) = Gi1(1)-2*D*ca(t(i))/incrX;
Gi1(end) = Gi1(end)+2*D*cb(t(i))/incrX;
Gi2 = g(x,t(i)+incrT); %Gi2 es el vector G en t(i)+k
Gi2(1) = Gi2(1)-2*D*ca(t(i)+incrT)/incrX;
Gi2(end) = Gi2(end)+2*D*cb(t(i)+incrT)/incrX;
K1 = -K*solheun(:,i)+Gi1;
K2 = -K*(solheun(:,i)+K1*incrT)+Gi2;
solheun(:,i+1) = solheun(:,i)+incrT/2*(K1+K2);
end
%% Representacion grafica comparativa
[Mt,Mx]= meshgrid(t,x);
% Error entre método del trapecio y Euler implicito
error1 = abs(soltrap-soleulerimp);
figure(1)
mesh(Mt,Mx,error1), shading flat;
title('Error entre los métodos del trapecio y Euler implícito')
xlabel('Tiempo (s)'); ylabel('Coordenada x (m)'); zlabel('Error')
colorbar
% Error entre método del trapecio y Euler explícito
error2 = abs(soltrap-soleuler);
figure(2)
mesh(Mt,Mx,error2), shading flat;
title('Error entre los métodos del trapecio y Euler explícito')
xlabel('Tiempo (s)'); ylabel('Coordenada x (m)'); zlabel('Error')
colorbar
% Error entre método del trapecio y Heun
error3 = abs(soltrap-solheun);
figure(3)
mesh(Mt,Mx,error3), shading flat;
title('Error entre los métodos del trapecio y Heun')
xlabel('Tiempo (s)'); ylabel('Coordenada x (m)'); zlabel('Error')
colorbar
% Error entre método del Euler implicito y Euler explicito
error4 = abs(soleulerimp-soleuler);
figure(4)
mesh(Mt,Mx,error4), shading flat;
title('Error entre los métodos de Euler implícito y Euler explícito')
xlabel('Tiempo (s)'); ylabel('Coordenada x (m)'); zlabel('Error')
colorbar
| Respecto a los métodos del trapecio y Euler implícito, se aprecia que difieren en pequeña magnitud en los primeros instantes de tiempo, pero según el tiempo va creciendo, los métodos convergen y su error disminuye a prácticamente 0. | 
|
| En cambio, en la comparativa del método del trapecio con el de Euler explícito no podemos observar el en los primeros instantes porque crece a medida que avanzamos en el tiempo, como se puede observar en el corte por el tiempo de 7 segundos, hasta valores muy altos, especialmente en el centro del tubo. Este orden de magnitud de 10^139 hace que no se observe en el gráfico el error en tiempos pequeños, que también existe. Esto se explica por la inestabilidad del método de Euler explícito que comentamos con anterioridad. | 
|
| La comparativa entre el método del trapecio y de Heun es muy similar a la anterior dado que el método de Heun también resultaba inestable. La unica diferencia está en el orden de magnitud que esta vez es de 10^226, mucho más alto que el método anterior. Además el error es más constante a lo largo del tubo que en la comparativa del trapecio con Euler explícito, es decir, hay menos diferencia entre el error en el punto medio y los extremos del tubo. | 
|
| Este último gráfico muestra la comparativa entre los dos métodos de Euler. El método implícito es el más exacto que hemos obtenido hasta ahora y se compara con el método explícito que resultó inestable. Por tanto se muestra una gráfica similar a las anteriores con un error de orden de magnitud 10^139 al igual que en la segunda comparación. La similitud entre la segunda comparación y la última corrobora el error entre el método del trapecio y de Euler implícito que obtuvimos al principio y que mostraba un error despreciable a medida que el tiempo avanzaba. | 
|
3 Conservación de la masa
Para demostrar el cumplimiento del principio de conservación de la masa a lo largo del tiempo, basta con integrar la ecuación diferencial de nuestro problema inicial en [math]x\in(0,L)[/math] con [math]L=8m[/math].
- [math]\frac{d}{dt}\displaystyle\int_{0}^{L}u(x,t)\,dx=\displaystyle\int_{0}^{L}u_{xx}(x,t)\,dx=u_x(x,t)|_0^L=0[/math]
La resolución de la segunda integral es inmediata pues equivale a la diferencia entre el valor de las condiciones frontera en [math]L=0[/math] y [math]L=8[/math], las cuales son nulas.
La primera integral representa la cantidad total de contaminante en el tubo. Por tanto, la expresión quedaría de la siguiente forma:
- [math]\frac{d}{dt}\displaystyle\int_{0}^{L}u(x,t)\,dx = \frac{d}{dt}(M_{total}) = 0[/math]
Lo que implica que [math]M_{total}[/math] será constante para [math]t\gt0[/math], con lo que queda demostrada la conservación de la masa.
Para obtener el valor de la cantidad total de contaminante, empleamos la resolución de la integral por el método del trapecio, el cual viene implementado en MATLAB con la función "trapz":
clear, close all
%% Datos del problema
a = 0; %Extremo derecho del tubo
b = 8; %Extremo izquierdo del tubo
t0 = 0;
tM = 12; %Tiempo estudiado
D = 1;
%% Funciones
% g : función del término independiente (x vector, t escalar)
g = @(x,t) x*0;
% ca : funcion f(t) de la condicion de contorno en x=a (t vector)
ca = @(t) t*0;
% cb : funcion f(t) de la condicion de contorno en x=b (t vector)
cb = @(t) t*0;
% u0 : forma analitica f(x) funcion valor inicial (x vector)
u0 = @(x) (2.*(x<=4))+(4.*(x>4));
%% Discretización del vector de espacio
incrX = 0.15;
N = round((b-a)/incrX);
x = linspace(a,b,N+1);
x = x';
%% Planteamiento del sistema a resolver U' = -K*U + G
% Matriz K
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2) = -2;
K(N+1,N) = -2;
K = K*(D/incrX^2);
% Condición inicial
U0 = u0(x);
%% Resolucion del sistema como PVI
% Discretización del vector de tiempos
incrT = incrX/4;
M = round((tM-t0)/incrT);
t = linspace(t0,tM,M+1);
%Inicializacion de la matriz de soluciones
sol = zeros(N+1,M+1);
sol(:,1)= U0;
% Método de Euler implícito
for i=1:M
Gi = g(x,t(i)); % Gi es el término independiente de la ecuación en el elemento n de la iteración (nulo en este caso)
Gi(1) = Gi(1)-2*D*ca(t(i))/incrX;
Gi(end) = Gi(end)+2*D*cb(t(i))/incrX;
Gi2 = g(x,t(i+1)); % Gi2 es el término independiente de la ecuación en el elemento n+1 de la iteración (nulo en este caso)
Gi2(1) = Gi2(1)-2*D*ca(t(i+1))/incrX;
Gi2(end) = Gi2(end)+2*D*cb(t(i+1))/incrX;
sol(:,i+1) = (eye(size(K))+incrT*K)\(sol(:,i)+incrT*Gi2);
end
%% Calculo de la primera integral de la expresión de la conservación de la masa
integral = zeros(1,length(t));
for i=1:length(t)
f=sol(:,i);
integral(i)=trapz(x,f);
end
%% Representacion grafica
plot(t,integral)
title('Evolución de la cantidad de contaminante con el tiempo')
xlabel('Tiempo (s)'); ylabel('Cantidad de contaminante (mol)');
disp(integral)
Con el código MATLAB, se obtiene la anterior gráfica en la que se representan la cantidad de contaminante en función del tiempo y se puede apreciar que el valor de la cantidad de contaminante permanece constante a lo largo del tiempo, tal y como se ha demostrado analíticamente en el comienzo del epígrafe.
4 Evolución de la concentración en el punto medio de la tubería
A continuación se analiza la evolución de la concentración de el punto medio de la tubería a través del siguiente código MATLAB. Como detalle a destacar, dado que el número de subintervalos para la aproximación numérica que se ha utilizado es 53 y, en consecuencia, el número de nodos a calcular será 54, para conocer la concentración en el punto medio tomaremos la media de las concentraciones de los nodos 27 y 28:
clear, close all
%% Datos del problema
a = 0; %Extremo derecho del tubo
b = 8; %Extremo izquierdo del tubo
t0 = 0;
tM = 12; %Tiempo estudiado
D = 1;
%% Funciones
% g : función del término independiente (x vector, t escalar)
g = @(x,t) x*0;
% ca : funcion f(t) de la condicion de contorno en x=a (t vector)
ca = @(t) t*0;
% cb : funcion f(t) de la condicion de contorno en x=b (t vector)
cb = @(t) t*0;
% u0 : forma analitica f(x) funcion valor inicial (x vector)
u0 = @(x) (2.*(x<=4))+(4.*(x>4));
%% Discretización del vector de espacio
incrX = 0.15;
N = round((b-a)/incrX);
x = linspace(a,b,N+1);
x = x';
%% Planteamiento del sistema a resolver U' = -K*U + G
% Matriz K
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2) = -2;
K(N+1,N) = -2;
K = K*(D/incrX^2);
% Condición inicial
U0 = u0(x);
%% Resolucion del sistema como PVI
% Discretización del vector de tiempos
incrT = incrX/4;
M = round((tM-t0)/incrT);
t = linspace(t0,tM,M+1);
%Inicializacion de la matriz de soluciones
sol = zeros(N+1,M+1);
sol(:,1)= U0;
% Método de Euler implícito
for i=1:M
Gi = g(x,t(i)); % Gi es el término independiente de la ecuación en el elemento n de la iteración (nulo en este caso)
Gi(1) = Gi(1)-2*D*ca(t(i))/incrX;
Gi(end) = Gi(end)+2*D*cb(t(i))/incrX;
Gi2 = g(x,t(i+1)); % Gi2 es el término independiente de la ecuación en el elemento n+1 de la iteración (nulo en este caso)
Gi2(1) = Gi2(1)-2*D*ca(t(i+1))/incrX;
Gi2(end) = Gi2(end)+2*D*cb(t(i+1))/incrX;
sol(:,i+1) = (eye(size(K))+incrT*K)\(sol(:,i)+incrT*Gi2);
end
%% Obtención de la concentración a lo largo del tiempo del punto medio
% N = 53 => nodos = 54 => Para el punto medio hacemos media de los elementos 27 y 28
Cpmedio = (sol(round(1+N/2),:)+sol(round(1+N/2+1),:))/2;
%% Representacion grafica
plot(t,Cpmedio)
title('Evolución de la concentración de contaminante en el punto medio con el tiempo')
xlabel('Tiempo (s)'); ylabel('Concentración (mol/(m^2*s))');
El gráfico que se obtiene con dicho código es el siguiente:
Al ser el último punto del tubo con concentración igual a 4, ya que los siguientes según el eje [math]x[/math] son puntos con concentración igual a 2, la concentración en este punto bajará rápidamente ya que es el primer punto que tiene que compensar la diferencia de concentraciones.
5 Estado estacionario
Para tiempos grandes, el contaminante se distribuye homogéneamente en el tubo a un valor de 3. Esto se puede apreciar en las gráficas de los métodos del trapecio y de Euler implícito. Las ecuaciones que debe satisfacer este estado estacionario son:
- [math] (P)\left\{\begin{matrix} u_{xx}=0 ,x\in\ (0,8), t\gt0\\u_x(0,t)=0, t\gt0\\u_x(8,t)=0, t\gt0\\u(x,0)=u_0, x\in\ (0,8)\end{matrix}\right. [/math]
La función de la solución será :[math]u{x}(x,t)=3[/math] para tiempos grandes.
Para ver la evolución de la concentración respecto al tiempo y la aproximación al estado estacionario, hemos calculado la diferencia entre la solución estacionaria y la solución en los puntos 0,1,2 y 10.
clc,clear all
%Datos
L=8; %Longitud de varilla
T=100; %Tiempo (usamos un valor alto, así nos aseguramos de alcanzar el estado estacionario)
D=1; %dato drl enunciado
%DISCRETIZACION TEMPORAL Y ESPACIAL
%Espacial
dx=0.15; % Paso en espacio
N=round(L/dx); %Número de subintervalos
x=linspace(0,L,N+1); % Vector de espacio
%Temporal
dt=dx/4; % Paso en tiempo
M=round(T/dt);%tiempo inicial t=0
t=linspace(0,T,M+1);
%t=0:dt:T; % Vector de tiempos
% Matriz de coeficientes. Aproximacion de -q*u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2;
K(N+1,N)=-2;
K=D*K/dx^2;
% Calculamos u0 condición inicial (datos del enunciado)
u0=[2*ones(1,round(N*4/8)),4*ones(1,round(N*4/8))]';
sol(1,:)=u0';
U=u0;
%Aplicamos método del trapecio
for j=1:length(t)-1
U=(eye(N+1)+(dt*K)/2)\((eye(N+1)-(K*dt)/2)*U);
sol(j+1,:)= U';
end
%una vez aplicada la iteración obtenemos la solucion estacionario a partir
%de la matriz de las soluciones
se=sol(length(t),:); %Solucion estacionaria
s0=sol(1,:); %Solución inicial (t=0)
pos1=round((1-0)/dt+1);
s1=sol(pos1,:);
pos2=round((2-0)/dt+1);
s2=sol(pos2,:);
pos10=round((10-0)/dt+1);
s10=sol(pos10,:);
%Ya hemos obtenido las soluciones de la matriz de soluciones, el siguiente
%paso es us análisis gráfico
hold on
plot(x,se,'-r');
plot(x,s0,'-c');
plot(x,s1,'-g');
plot(x,s2,'-y');
plot(x,s10,'-b');
legend('u(Estacionaria)','u(Inicial (t=0))','u(1 segundo)','u(2 segundos)','u(10 segundos)');
xlabel('Posición en la varilla')
ylabel('Concentracion')
hold off
Se observa cómo la solución se aproxima al estado estacionario con el tiempo.
5.1 Tiempo en alcanzar la solución estacionaria
El estado estacionario es el valor de concentración u = 3 como podemos comprobar en el gráfico y obteniendo los valores en cualquier punto de x para un valor de tiempo alto. Para encontrar el tiempo en que se alcanza este valor consideramos un error del 5%, que corresponderá a 0,15. Por tanto buscaremos el tiempo a partir del cual todos los puntos de la barra tienen su concentración comprendida en el intervalo [2,85-3,15]. Ayudándonos por la matriz solución y sabiendo que los extremos de la barra corresponden a los valores máximo y mínimo de la barra en todo instante (ver gráfico), tendiendo ambos asintóticamente al 3 desde los valores 2 y 4, basta con comprobar cuando los valores en los extremos llegan a entrar en el intervalo. Dado que 3 es el valor medio entre los valores iniciales en ambos extremos 2 y 4, la aproximación se produce de forma simétrica y en el punto de columna 364 de la matriz solución encontramos los valores dentro del intervalo pedido. Pidiendo a matlab el tiempo en que se produce este elemento con el comando t(364) obtenemos el valor de tiempo 13.77233 segundos, que es el tiempo que tarda en alcanzar el estado estacionario con error del 5%.
La solución cambia ligeramente si ∆x se divide entre diez ya que la aproximación al resultado real será mejor. Pero la variación es muy pequeña dado que ya teníamos un ∆x suficientemente pequeño para la resolución. Lo encontramos en la columna 3695 y nos da un valor de 13.8523 segundos.
