Filtration modelisation pdepe derivation problem
    10 views (last 30 days)
  
       Show older comments
    
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 . 
0 Comments
Answers (1)
  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.
0 Comments
See Also
Categories
				Find more on Assembly 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!
