ANEXO II
Scripts-Solución Analítica
Código de Matlab para el primer escenario (fuente radioactiva al nivel del suelo) de
la solución analítica
%Caso 12
%Dispersión de pluma/Estado Estacionario en dirección x a una
velocidad
%constante del viento
%C(x,y,z)=exp^(-0.5*(y^2/thetay^2
+z^2/thetaz^2)*Qm/pi*thetaz*thetay*u
clc
clear all
close all
CASE = input('Ingrese el escenario: ','s');
%Para cada par de y existe un z
Qm=1000; %Flujo másico [kg/s]
u=20; %Velocidad del viento [m/s]
%Se toma en cuenta los vectores
[y,x]=meshgrid(1e-4:10:2000,1e-4:10:2000);
z =50;
%Para la condición ambiental:
sigma_y = fsigma_y(CASE,x); %[m]
sigma_z = fsigma_z(CASE,x); %[m]
exponencial = exp(-0.5*(((y.^2)./(sigma_y.^2))
+((z^2)./(sigma_z.^2))));
C =exponencial.*(Qm./(pi*sigma_z.*sigma_y.*u));
figure(1)
surf(x,y,C)
figure(2)
contourf(x,y,C,100,'LineColor','none')
xlabel('Distancia [m]')
ylabel('Distancia [m]')
zlabel('Concentración [kg/m^3]')
hold on
contourf(x,-y,C,100,'LineColor','none')
axis([0 1990 -200 200])
grid on
59
Código de Matlab para el segundo escenario (fuente radioactiva a una determinada
altura “H”) de la solución analítica
%Caso 13
%Dispersión de pluma/Estado Estacionario en dirección x a una
velocidad
%constante del viento y determinada altura
%C(x,y,z)=(exp(-0.5*(((y.^2)./(sigma_y.^2))))*exp(-0.5*(((z-
H)^2)./(sigma_z.^2)))+exp(-
0.5*(((z+H)^2)./(sigma_z.^2))))*Qm/pi*thetaz*thetay*u
clc
clear all
close all
CASE = input('Ingrese el escenario: ','s');
%Para cada par de y existe un z
Qm=1000; %Flujo másico [kg/s]
u=50; %Velocidad del viento [m/s]
%Se toma en cuenta los vectores
[y,x]=meshgrid(1e-4:10:2000,1e-4:10:2000);
z =50; %Altura de la chimenea
H= 10; %Altura sobre la que está la chimenea
%Para la condición ambiental:
sigma_y = fsigma_y(CASE,x); %[m]
sigma_z = fsigma_z(CASE,x); %[m]
exponencial1 = exp(-0.5*(((y.^2)./(sigma_y.^2))));
exponencial2 = exp(-0.5*(((z-H)^2)./(sigma_z.^2)))+exp(-
0.5*(((z+H)^2)./(sigma_z.^2)));
C =exponencial1.*exponencial2.*(Qm./(pi*sigma_z.*sigma_y.*u));
figure(1)
surf(x,y,C)
figure(2)
contourf(x,y,C,100,'LineColor','none')
xlabel('Distancia [m]')
ylabel('Distancia [m]')
zlabel('Concentración [kg/m^3]')
hold on
contourf(x,-y,C,100,'LineColor','none')
axis([0 1990 -200 200])
grid on
60
Código para la determinación de coeficiente de dispersión "𝝈𝒚 " de la solución
analítica
function [sigma_y]=fsigma_y(CASE,x)
%CASE = input('Ingrese el caso','s');
ks1 = {'RA','RB','RC','RD','RE','RF'};
v1 = {0.22,0.16,0.11,0.08, 0.06, 0.04};
ks2 = {'UA','UB','UC','UD','UE','UF'};
v2 = {0.32,0.32,0.22,0.16, 0.11, 0.11};
M1 = containers.Map(ks1,v1);
M2 = containers.Map(ks2,v2);
if CASE(1)=='R'
c = M1(CASE);
sigma_y=c*x.*(1+0.0001*x).^-0.5;
elseif CASE(1)=='U'
c = M2(CASE);
sigma_y=c*x.*(1+0.0004*x).^-0.5;
else
disp('Este caso no existe')
end
end
61
Código para la determinación de coeficiente de dispersión "𝝈𝒛 " de la solución
analítica
function [sigma_z]=fsigma_z(CASE,x)
%CASE = input('Ingrese el caso','s');
ks1 = {'RA','RB','RC','RD','RE','RF'};
v1 = {[0.20,0],[0.12,0],[0.08, 0.0002],[0.06,0.0015],[0.03,
0.0003],[0.016, 0.0003]};
ks2 = {'UA','UB','UC','UD','UE','UF'};
v2 = {[0.24,0.0001],[0.24,0.0001],[0.20,0],[0.14,0.0003],[0.08,
0.0015], [0.08, 0.0015]};
M1 = containers.Map(ks1,v1);
M2 = containers.Map(ks2,v2);
if CASE(1)=='R'
c = M1(CASE);
if CASE=='RE' | CASE=='RF'
sigma_z=c(1)*x.*(1+c(2)*x).^-1;
else
sigma_z=c(1)*x.*(1+c(2)*x).^-0.5;
end
elseif CASE(1)=='U'
c = M2(CASE);
if CASE=='UA' | CASE=='UB'
sigma_z=c(1)*x.*(1+c(2)*x).^0.5;
else
sigma_z=c(1)*x.*(1+c(2)*x).^-0.5;
end
else
disp('Este caso no existe')
end
end
62
ANEXO III
Scripts-Solución Numérica
Código de Matlab de la solución numérica por métodos de diferencias finitas y Euler
clc
clear all
close all
%Método de Diferencias Finitas y Euler
%Definición de constantes
Qm = 10; % Flujo de liberación de la fuente [kg/s]
k = 0.086; % Coeficiente de difusión [m^2/s]
u = 1; % Velocidad del aire [m/s]
tf=1; % Tiempo final[s]
%Criterio de estabilidad
K=0.9; % Coeficiente de proporcionalidad
N = 100; % Número de nodos
L = 0.8; % Longitud [m]
dx = L/(N-1);
dt=K*(dx^2)/(2*k);
Co = 2; % Concentración inicial [kg/m^3]
%Definición de constantes
A = u*dt/(2*dx) + k*dt/(dx^2);
B = 1 - 2*k*dt/(dx^2);
D = -u*dt/(2*dx)+k*dt/(dx^2);
C = zeros(N,N);
x = linspace(0,L,N);
t = linspace(0,tf,N);
% i se refiere al tiempo
% j se refiere a la posición
C(1,:)=Co;
for i=1:N-1
for j=2:N-1
C(i+1,j)=A*C(i,j-1)+B*C(i,j)+D*C(i,j+1);
end
end
dist = [2,round(N/2),N-1];
titulo = sprintf('Solución numérica (%6.0f puntos de malla)',N);
subplot(2,2,1)
plot(t,C(:,dist(1)),'LineWidth',1.2)
title(titulo)
63
xlabel('Tiempo [s]')
ylabel('Concentración [kg/m^3]')
leyenda1 = 'Distancia %2.4f m';
leyenda1 = sprintf(leyenda1,x(dist(1)));
legend(leyenda1)
grid on
hold on
subplot(2,2,2)
plot(t,C(:,dist(2)),'LineWidth',1.2)
title(titulo)
xlabel('Tiempo [s]')
ylabel('Concentración [kg/m^3]')
leyenda1 = 'Distancia %2.4f m';
leyenda1 = sprintf(leyenda1,x(dist(2)));
legend(leyenda1)
grid on
subplot(2,2,3)
plot(t,C(:,dist(3)),'LineWidth',1.2)
title(titulo)
xlabel('Tiempo [s]')
ylabel('Concentración [kg/m^3]')
leyenda1 = 'Distancia %2.4f m';
leyenda1 = sprintf(leyenda1,x(dist(3)));
legend(leyenda1)
grid on
hold off
tiem = [2,round(N/5), round(N*2/5),round(N*3/5), round(N*4/5), N-
1];
figure(2)
legend2 = [];
for i = tiem
plot(x,C(i,:),'LineWidth',1.2)
hold on
legend2=[legend2; sprintf('Tiempo %2.4f seg',t(i)) ];
end
title(titulo)
xlabel('Distancia [m]')
ylabel('Concentración [kg/m^3]')
legend(legend2)
grid on
64
Código de Matlab de la solución numérica por el método de derivadas parciales
clc
clear all
close all
%Método Función integrada de Matlab
C.k = 0.086; %Constante de Difusion
C.uv = 1; %Velocidad del viento
C.L = 1; %Longitud
C.Co = 2; %Concentración inicial
C.tf=1; %Tiempo final
N=100;
x = linspace(0,C.L,N); %Número de mallado
t = linspace(0,C.tf,N);
m = 0;
eqn = @(x,t,Con,dCondx) ecuacionPDE(x,t,Con,dCondx,C); %EDP
ic = @(x) condicioninicial(x,C); %Ecuación condición
inicial
Con = pdepe(m,eqn,ic,@condicionfrontera,x,t); %Ecuación
condicones de frontera
%pdepe (ecuación
integradaMatlab)
dist = [2,round(N/2),N-1];
titulo = sprintf('Solución numérica (%6.0f puntos de malla)',N);
subplot(2,2,1)
plot(t,Con(:,dist(1)),'LineWidth',1.2)
title(titulo)
xlabel('Tiempo [s]')
ylabel('Concentración [kg/m^3]')
leyenda1 = 'Distancia %2.4f m';
leyenda1 = sprintf(leyenda1,x(dist(1)));
legend(leyenda1)
grid on
hold on
subplot(2,2,2)
plot(t,Con(:,dist(2)),'LineWidth',1.2)
title(titulo)
xlabel('Tiempo [s]')
ylabel('Concentración [kg/m^3]')
leyenda1 = 'Distancia %2.4f m';
leyenda1 = sprintf(leyenda1,x(dist(2)));
legend(leyenda1)
grid on
subplot(2,2,3)
plot(t,Con(:,dist(3)),'LineWidth',1.2)
title(titulo)
xlabel('Tiempo [s]')
65
ylabel('Concentración [kg/m^3]')
leyenda1 = 'Distancia %2.4f m';
leyenda1 = sprintf(leyenda1,x(dist(3)));
legend(leyenda1)
grid on
hold off
tiem = [2,round(N/5), round(N*2/5),round(N*3/5), round(N*4/5), N-
1];
figure(2)
legend2 = [];
for i = tiem
plot(x,Con(i,:),'LineWidth',1.2)
hold on
legend2=[legend2; sprintf('Tiempo %2.4f seg',t(i)) ];
end
title(titulo)
xlabel('Distancia [m]')
ylabel('Concentración [kg/m^3]')
legend(legend2)
grid on
function [c,f,s] = ecuacionPDE(x,t,Con,dCondx,C)
k = C.k;
uv = C.uv;
c = 1; %Derivada temporal
f = k*dCondx; %Segunda derivada espacial
s = -(uv)*dCondx; %Primera derivada espacial
end
function u0 = condicioninicial(x,C)
u0 = C.Co;
end
function [pl,ql,pr,qr] = condicionfrontera(xl,Conl,xr,Conr,t)
pl = Conl; %Condición de borde izquierda
ql = 0; %Valor de la condición de borde izquierda
pr = Conr; %Condición de borde derecha
qr = 0; %Valor de la condición de borde derecha
end
66
