Filtration modelisation pdepe derivation problem

1 Ansicht (letzte 30 Tage)
Juliette Parrott
Juliette Parrott am 12 Jan. 2021
Beantwortet: Josh Meyer am 27 Jan. 2021
Hello ,
I would like to modelize a filtration using Pdepe. I have already done a part of it : i have my main function f and I added function D in it so D is now counted in F. I would like to add an other parametre : J. J should be function of phi and ruled by this equation : J*phi-d(phi)*dphidx=0
How can I include this function ? I have tried with "diff" but it is not working. Mu main problem is on the dphidx part : is there a method to calculate this part?
here is my code with a J constant :
clear
close all
%initiation
J=110/3600; % L.s^-1.m^-2
L=7*10^-3; %p115
Ca=0.02; % on choisi la betonite - page 86 -(VFA*Ca=1,8 et VFA = 85)
N=100;
I=100 ;
m=0;
phi=linspace(0,0.7,100);
xmesh=linspace(0,0.5e-6,50);
tspan=linspace(0,5e-4,10);
for i=1:100
D(i)=fD(phi(i))
end
figure(1)
plot(phi,D)
axis([0 0.7 0 5e-9]);
title('évolution du coefficient de diffusion')
xlabel('fraction volumique')
ylabel('coefficient de diffusion')
%stop
sol=pdepe(m,@pdefun,@icfun,@bcfun,xmesh,tspan);
figure(2)
surf(xmesh,tspan,sol)
title('evolution du profil de concentration en fonction du temps et épaisseur')
xlabel('épaisseur')
zlabel('concentration')
ylabel('temps')
%fonction function D=fD(phi,phic)
function D=fD(phi)
%phi=phi/10;
a=11.8;
b=2.9;
g=-4.7e2;
d=-2.1e3;
e=-44.6e3;
f=26.3e3;
phic=0.5723;
Rp=90*10^9;
mus=10^-3;
Vp=0.00305;
if phi<phic
den=(1/(a*phi+b))+(1/(g*phi+d))+(1/(e*phi+f));
po=exp(1/den);
numpophi=(a/((a*phi+b)^2))+(g/((g*phi+d)^2))+(e/((e*phi+f)^2));
dpodphi=po*(numpophi/(den)^2);
h=(6+4*phi^(5/3))/(6-9*phi^(1/3)+9*phi^(5/3)-6*phi^2);
D=(1/(6*pi*mus*Rp*h))*Vp*dpodphi;
else
D=(phi-phic)*5e-9/2e-2;
%D=5e-9;
end
end
function C0=icfun(x)
Ca=0.02;
C0=Ca;
end
function [pl,ql,pr,qr]=bcfun(xl,Cl,xr,Cr,t);
Ca=0.02;
pl=0;
ql=1;
pr=Cr-Ca;
qr=0;
end
function [c,f,s]=pdefun(x,t,C,dCdx)
J=110/3600;
s=0;
c=1;
f=J*C+fD(C)*dCdx;
end
Thanks in advance for your help .

Antworten (1)

Josh Meyer
Josh Meyer am 27 Jan. 2021
Is phi supposed to be a function of x? You currently have
phi=linspace(0,0.7,100);
and
xmesh=linspace(0,0.5e-6,50);
So, there appears to be no relation between the two, and the vectors have different lengths. However, if you calculated phi as a function of x, for example
phi=xmesh.^2;
Then you would be able to use finite differences to approximate d(phi)/dx on the grid:
dphidx = diff(phi)./diff(xmesh);
Note that if phi is not a function of x, then d(phi)/dx = 0.

Kategorien

Mehr zu Partial Differential Equation Toolbox finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by