Borrar filtros
Borrar filtros

Why my graph not same as research paper?

2 visualizaciones (últimos 30 días)
nur
nur el 8 de Jun. de 2022
Comentada: Asif Solanki el 25 de Ag. de 2022
Hi, i want to ask why i did not get same graph as above?
This is my code that i have try:
%for t1=5%
x=0:0.2:10;
m=0.1;
C_0=1.0;
u_0=0.25;
D_0=0.45;
w_0=0.001;
p_0=0.02;
t1=5;
a=((((u_0)^2)/(4*D_0))+p_0);
b=((u_0)^2)/(4*D_0);
d=(1/(sqrt(2)*m));
e1=(log(2+sqrt(2)+(4+3*(sqrt(2)))*(tanh(m*t1/2))));
e2=(log(2+sqrt(2)-sqrt(2)*tanh(m*t1/2)));
T=d*(e1-e2);
f=((x-((((u_0)^2)+4*p_0*D_0)^1/2)*T)/(2*sqrt((D_0)*T)));
g=((x+((((u_0)^2)+4*p_0*D_0)^1/2)*T)/(2*sqrt((D_0)*T)));
h=((((u_0)-((((u_0)^2)+4*p_0*D_0)^1/2))*x)/(2*D_0));
i=((((u_0)+((((u_0)^2)+4*p_0*D_0)^1/2))*x)/(2*D_0));
j=((x-(u_0)*T)/(2*sqrt((D_0)*T)));
k=((x+(u_0)*T)/(2*sqrt((D_0)*T)));
l=(((u_0)*x)/(D_0));
n=exp(h);
o=erfc(f);
q=exp(i);
r=erfc(g);
s=exp(-p_0*T);
v=erfc(j);
y=erfc(k);
z=exp(l);
F=(((1/2).*n.*o)+((1/2).*q.*r));
G=(s.*(1-(1/2).*v-(1/2).*z.*y));
C1=(((w_0)/(p_0))+(C_0-((w_0)/(p_0))))*F-(((w_0)/(p_0))*G);
plot(x,C1)
hold on
%for t2=10%
x=0:0.2:10;
m=0.1;
C_0=1.0;
u_0=0.25;
D_0=0.45;
w_0=0.001;
p_0=0.02;
t2=10;
a=(((u_0)^2)/(4*D_0)+p_0);
b=((u_0)^2)/(4*D_0);
d=(1/(sqrt(2)*m));
e1=(log(2+sqrt(2)+(4+3*(sqrt(2)))*(tanh(m*t2/2))));
e2=(log(2+sqrt(2)-sqrt(2)*tanh(m*t2/2)));
T=d*(e1-e2);
f=((x-((((u_0)^2)+4*p_0*D_0)^1/2)*T)/(2*sqrt((D_0)*T)));
g=((x+((((u_0)^2)+4*p_0*D_0)^1/2)*T)/(2*sqrt((D_0)*T)));
h=((((u_0)-((((u_0)^2)+4*p_0*D_0)^1/2))*x)/(2*D_0));
i=((((u_0)+((((u_0)^2)+4*p_0*D_0)^1/2))*x)/(2*D_0));
j=((x-(u_0)*T)/(2*sqrt((D_0)*T)));
k=((x+(u_0)*T)/(2*sqrt((D_0)*T)));
l=(((u_0)*x)/(D_0));
n=exp(h);
o=erfc(f);
q=exp(i);
r=erfc(g);
s=exp(-p_0*T);
v=erfc(j);
y=erfc(k);
z=exp(l);
F=(((1/2).*n.*o)+((1/2).*q.*r));
G=(s.*(1-(1/2).*v-(1/2).*z.*y));
C2=(((w_0)/(p_0))+(C_0-((w_0)/(p_0))))*F-(((w_0)/(p_0))*G);
plot(x,C2)
hold off
hold on
%for t3=15%
x=0:10;
m=0.1;
C_0=1.0;
u_0=0.25;
D_0=0.45;
w_0=0.001;
p_0=0.02;
t3=15;
a=(((u_0)^2)/(4*D_0)+p_0);
b=((u_0)^2)/(4*D_0);
d=(1/(sqrt(2)*m));
e1=(log(2+sqrt(2)+(4+3*(sqrt(2)))*(tanh(m*t3/2))));
e2=(log(2+sqrt(2)-sqrt(2)*tanh(m*t3/2)));
T=d*(e1-e2);
f=((x-((((u_0)^2)+4*p_0*D_0)^1/2)*T)/(2*sqrt((D_0)*T)));
g=((x+((((u_0)^2)+4*p_0*D_0)^1/2)*T)/(2*sqrt((D_0)*T)));
h=((((u_0)-((((u_0)^2)+4*p_0*D_0)^1/2))*x)/(2*D_0));
i=((((u_0)+((((u_0)^2)+4*p_0*D_0)^1/2))*x)/(2*D_0));
j=((x-(u_0)*T)/(2*sqrt((D_0)*T)));
k=((x+(u_0)*T)/(2*sqrt((D_0)*T)));
l=(((u_0)*x)/(D_0));
n=exp(h);
o=erfc(f);
q=exp(i);
r=erfc(g);
s=exp(-p_0*T);
v=erfc(j);
y=erfc(k);
z=exp(l);
F=(((1/2).*n.*o)+((1/2).*q.*r));
G=(s.*(1-(1/2).*v-(1/2).*z.*y));
C3=(((w_0)/(p_0))+(C_0-((w_0)/(p_0))))*F-(((w_0)/(p_0))*G);
plot(x,C3)
hold on
legend({'t=5','t=10','t=15'},'Location','southwest')
xlabel('Distance,x (meter)')
ylabel('Concentration,C(x,t)')
%ylim([0, 1])
xlim([0, 10])
  6 comentarios
Torsten
Torsten el 8 de Jun. de 2022
The only thing that helps:
Compare the formulas from the article and yours.
Torsten
Torsten el 8 de Jun. de 2022
I corrected the three errors in your code detected by @SALAH ALRABEEI
Looks better now, but concentrations still become negative. So your search must go on.

Iniciar sesión para comentar.

Respuesta aceptada

Torsten
Torsten el 8 de Jun. de 2022
Editada: Torsten el 8 de Jun. de 2022
Sometimes it's better to start anew:
x=(0:0.1:10).';
C0 = 1.0;
tt = [5.0,10.0,15,0];
u0 = 0.25;
D0 = 0.45;
m = 0.1;
gamma0 = 0.02;
mu0 = 0.001;
C = zeros(numel(x),numel(tt));
for i=1:numel(tt)
t = tt(i);
T = 1/(sqrt(2)*m)*(log(2+sqrt(2)+(4+3*sqrt(2))*tanh(m*t/2))-...
log(2+sqrt(2)-sqrt(2)*tanh(m*t/2)));
F = 0.5*exp((u0-sqrt(u0^2+4*gamma0*D0))*x/(2*D0)).*...
erfc((x-sqrt(u0^2+4*gamma0*D0)*T)./(2*sqrt(D0*T)))+...
0.5*exp((u0+sqrt(u0^2+4*gamma0*D0))*x/(2*D0)).*...
erfc((x+sqrt(u0^2+4*gamma0*D0)*T)./(2*sqrt(D0*T)));
G = exp(-gamma0*T)*(1-0.5*erfc((x-u0*T)./(2*sqrt(D0*T)))-...
0.5*exp(u0*x/D0).*erfc((x+u0*T)./(2*sqrt(D0*T))));
C(:,i) = mu0/gamma0 + (C0-mu0/gamma0)*F - mu0/gamma0*G;
end
plot(x,C)
  2 comentarios
nur
nur el 9 de Jun. de 2022
thank you so much!!! this graph exactly same as in the research paper. :)
Asif Solanki
Asif Solanki el 25 de Ag. de 2022
Dear Sir,
My name is Asif Solanki I am a student persuing master's program in Italy. I am currently working on a small project where I do need to find the concentration of the pollutants using the Bear analytical method.
Therefore, would you like to assist me to find the solution?
Thank you very much
Best regards
Asif Solanki

Iniciar sesión para comentar.

Más respuestas (2)

Sam Chak
Sam Chak el 8 de Jun. de 2022
Hi @nur
One of the ways to find out is to determine the equilibrium points from the advection-dispersion equation.
I haven't checked the long equations. Can you verify if you plotted C or ? Sometimes, the transformations can be a little tricky.
  2 comentarios
nur
nur el 8 de Jun. de 2022
I plotted C. I hope you can help me ;(
Sam Chak
Sam Chak el 8 de Jun. de 2022
Editada: Sam Chak el 8 de Jun. de 2022
Hi @nur
Try to make a little bit of value-added effort to the problem.
Another way is to plot the concentration C from the numerical solution of the advection-dispersion differential equation.
This way you can compare with the analytical solution obtained from the paper. Bear in mind that sometimes misprints can occur due to authors' mistake, or the production crew's mistake. So, what was shown on the paper might not be truly the analytical solution. Therefore, yes you have to check.
But I'd suggest you to take numerical solution approach because that is directly from the governing advection-dispersion law. The analytical solution, which is "human-processed equation" from the law.
Edit: Hey, check this out:
One of the objectives is to Solve the one-dimensional advection dispersion equation using ode45 or ode15s.

Iniciar sesión para comentar.


SALAH ALRABEEI
SALAH ALRABEEI el 8 de Jun. de 2022
You have two mistakes here
s=(exp(p_0)*T);
The three s in the code should be this
s=(exp(-p_0*T));
  3 comentarios
nur
nur el 8 de Jun. de 2022
thank youuu it much better but can i know why the end line of graph not end at 0 on y axis?
SALAH ALRABEEI
SALAH ALRABEEI el 8 de Jun. de 2022
The reason is that the curves have data less than zero. Unlike those in the paper.
Anyway, you still can do this by adding this
axis([0 10 0 1])

Iniciar sesión para comentar.

Categorías

Más información sobre Numerical Integration and Differential Equations en Help Center y File Exchange.

Etiquetas

Productos

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by