i want the area under a curve which varies in time.

1 次查看（过去 30 天）
Alfonso
Alfonso 2023-1-18
评论： Torsten 2023-1-20
The code is:
function ML_Diff
clc
clear all
%This is the fucking file of the exam exercise
%Explicit Euler Method is used for solving the following
%initial-boundary-value problem:
%Ut+v1*Ux=alpha1Uxx Ut+v2*Ux=alpha2Uxx
%U(z,0)=0 %Cond iniziale
%U(0,t)=g1(t) U(L,t)=g2(t) %Cond. Al bordo
%Problem parameters
nt = 100; %numero punti griglia asse x
nz = 150; %numero punti griglia asse y
tc = 1; %è il tempo caratteristico
T = 1000; %è il tempo adimensionalizzato t/tc
Z = 50; %profondità di tutto il multistrato
dZ = Z/nz; %è l'intervallo spaziale tra 2 punti zi
ni=1*10^-6;
%---------------------layer 1----------------------
z1 = Z/2; %z1 è la z barra del layer 1 ovvero z1=(L/2)/L
alpha1= 0.05;%0.05
v1 = 0.050;%0.050;
Pe1 = (v1*z1)/alpha1;
Dm1=(z1^2)/Pe1;
Re1=(v1*z1)/ni;
%---------------------Interface--------------------
zh=Z/2; %profondità dell'interfaccia.
%---------------------layer 2----------------------
z2 = Z/2; %z1 è la z barra del layer 1 ovvero z1=(L/2)/L
alpha2=0.1;
v2 = 0.050;%0.005;
Pe2 = (v2*z2)/alpha2;
Dm2=(z2^2)/Pe2;
Re2=(v2*z2)/ni;
%Check of problem parameters
if any([alpha1 alpha2 nz-2 nt-2 z1 z2]<=0)
error('Check problem parameters')
end
%Stability
st = 1;
while (2*max(alpha1,alpha2)*(T/nt/dZ^2)+(T/nt/dZ)*abs(max(v1,v2)))>st
nt=nt+1;
end
%Post-stabilità andiamo allora ad esplicitare i seguenti termini:
dT=T/nt; % Non abbiamo inserito anche il delta_Z perchè questo non
%cambia, cioè è lo stesso definito ad inizio codice.
r=(dT/dZ^2); s=(dT/dZ);
%Inizialization
h=(zh/dZ)+1; h=double(uint16(abs(h))); %index of interface
z=linspace(0,Z,nz+1); t=linspace(0,T,nt+1);
zh=z(h);
u=getPisciazz(z,nz,Z); U=u; %Fissiamo la nostra funzione concentrazione(u) pari prorpio alla
p1=alpha1*r;
d1=v1*s;
p2=alpha2*r;
d2=v2*s;
%Subfunction [Per le condizioni iniziali]
function Pisciazz = getPisciazz(z,nz,Z)
Pisciazz(1:nz+1,1)=0;
%Pisciazz(1:h,1)=0 %Cond iniziale asse z, 1°Strato
%Pisciazz(h+1:nz+1,1)=0 %Cond iniziale asse z, 2°Strato
end
%Subfunction [Per bordo z=0]
function VarTempo = getVarTempo(t)
%VarTempo=2
Dt=100; C0=50;
if t<=Dt
VarTempo= C0;
else
VarTempo= 0;
end
end
%Subfunction [Per bordo z=L]
function NoemannNothing = getNoemannNothing(t)
NoemannNothing= u(end-1);
%u(end-u)
%noemann nulla du/dz=0 u(end)=u(end-u)
%end
%Finite Difference Method [Forward for the time and Central for the space]
for j=2:nt+1
u(2:h-1)=u(2:h-1)-d1*(u(3:h)-u(2:h-1))+p1*(u(3:h)-2*u(2:h-1)+u(1:h-2));
u(h+1:end-1)=u(h+1:end-1)-d2*(u(h+2:end)-u(h+1:end-1))+p2*(u(h+2:end)-2*u(h+1:end-1)+u(h:end-2));
u(h)=(Dm1*(-4*u(h-1)+u(h-2))-(Dm2*(4*u(h+1)-u(h+2))))/(v1*2*dZ-v2*2*dZ-3*Dm1-3*Dm2);
%u(h)=(((4*u(h-1)-u(h-2)))*-Dm1-(4*u(h+1)-u(h+2))*Dm2)/((v1-v2)*2*dZ-3*(Dm1+Dm2));
u(1)=getVarTempo(t(j));
u(end)=getNoemannNothing(t(j));
plot(u,z,'r*:',U,z,'k');
set(gca,'YDir','reverse');
xlabel('C(z,t)'); ylabel('z');
title(['t=',num2str(t(j)),'seconds']);
axis([0 35 0 50]);
pause(0.01);
end
legend(['t=',num2str(T)],'t=0');
end
Please help me faster! I have to send the project to my teacher!
  13 个评论
显示 11更早的评论
Alfonso
Alfonso 2023-1-20
The fact Is that the area under the curve Is Just the Mass of contaminant, i've ran the code and discovered that the area Is not the same during che process above, but Is should be.
I think It Is a provlem for the diffusivity approximation terme, so maybe its wrong the centrale approximation, dont you think Is Is?
Torsten
Torsten 2023-1-20
Didn't you read my answer ? The mass of the contaminant changes with time (increases) since you inject more and more of it over a period of 100 s. So how should it be the same during the process ?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Electromagnetics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by