Filtration modelisation pdepe derivation problem

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 .

Answer 1

[Josh Meyer]

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.