Diferencia entre revisiones de «Difusión de una sustancia contaminante (Grupo 6)»
(→Resolución por el método de Euler explícito) |
|||
| Línea 234: | Línea 234: | ||
[[Archivo:Grafico_euler1b.png|1200px|center]] | [[Archivo:Grafico_euler1b.png|1200px|center]] | ||
| − | 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 * | + | 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). |
| + | |||
==== Resolución por el método de Euler implícito ==== | ==== Resolución por el método de Euler implícito ==== | ||
El código MATLAB empleado para su resolución es el siguiente: | El código MATLAB empleado para su resolución es el siguiente: | ||
Revisión del 15:33 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 de la barra, 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 del 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))+incrX/2*K)\(sol(:,i)+incrX/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:
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)+incrX*(-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))+incrX*K)\(sol(:,i)+incrX*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:
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:
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))+incrX/2*K)\(soltrap(:,i)+incrX/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)+incrX*(-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))+incrX*K)\(soleulerimp(:,i)+incrX*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
| 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. | 
|
| Respecto a los métodos del trapecio y Euler explícito, al ser el método explícito por el tamaño de paso tomado e 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. | 
|
| 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. | 
|
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))+incrX*K)\(sol(:,i)+incrX*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))+incrX*K)\(sol(:,i)+incrX*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 de la barra 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.
