Filtration modelisation pdepe derivation problem

1 view (last 30 days)
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 .

Answers (1)

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

Categories

Find more on Partial Differential Equation Toolbox in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!