# Filtration modelisation pdepe derivation problem

1 view (last 30 days)
Juliette Parrott on 12 Jan 2021
Answered: Josh Meyer on 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

Josh Meyer on 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.