pdepe help! The solution gives 0...

Hello,

My goal is to find the length L of the cylindrical column by solving this eqn:

Initial condition and Boundary conditions are as followed:

I attempted to find L by first solving the problem to find the eqn of C(t,z). Then by using (given value: Cf). I can find L. However, when I ran this code (as shown below), I get 0..... I have no idea where my code has gone wrong. If anyone could help me with this, I'll appreciate a lot.

Here is the code I have:

% Listing variables 
Cin=3.5;%g/L 
q_max=15;%g/L 
V=480; %L 
yield=91; % 
k_eq=7.7*10^(-4); % mol/m^3 
Cf=yield*Cin; 
epsilon=0.36; 
u0=0.00194; %m/s 
d_p=85*10^-9; % m 
Rho=1000; %g/L 
eta=1; % mPa s

% Find the missing variables 
Da=findDa(d_p,u0,epsilon,Rho,eta); %function 
constant=findConstant(epsilon, q_max, k_eq); 
% L=findL(u, Da,Cf); ----> goal

x = linspace(0,5,20); 
t = linspace(0,5,10);

% Solve 
m = 1; %cylindrical 
eqn = @(x,t,u,dudx) setPDE(x,t,u,dudx,Da,constant,u0); 
IC = @(x) findIC(x); 
bcfcn = @(xl,ul,xr,ur,t) setBC(xl,ul,xr,ur,t,Cin); 
sol = pdepe(m,eqn,IC,bcfcn,x,t);
% Extract the soln? 
u1 = sol(:,:,1);

surf(x,t,u) 
title('Numerical solution') 
xlabel('Distance z') 
ylabel('Time t') 
figure 
plot(x,u(end,:)) 
title('Solution at t = ') 
xlabel('Distance z') 
ylabel('C(z, )')

%-----------------Find variables-------------------% 
function constantVal=findConstant(epsilon, q_max, k_eq) 
    constantVal=(1+(1-epsilon)/epsilon*q_max*k_eq); 
end 
function Daval=findDa(d_p,u0,epsilon,Rho,eta) 
    Re=u0*d_p*Rho/eta; 
    Daval=d_p*u0*epsilon/(0.339+0.033*Re^0.48); 
end

%--------------------PDE---------------------------% 
function [c,f,s]=setPDE(x,t,u,dudx,Da,constant,u0) 
    c=1; 
    f=(Da/constant)*dudx; 
    s=(-u0/constant)*dudx; 
end 
%--------------------IC----------------------------% 
function C0=findIC(x) 
    C0=0; 
end 
%--------------------BC----------------------------% 
function [pl,ql,pr,qr]=setBC(xl,ul,xr,ur,t,Cin) 
    pl = ul-Cin; 
    ql = 0; 
    pr = 0; 
    qr = 1; 
end

Answer 1

[Torsten]

IC = @(x) findIC(x,Cin); 

...
    
function C0=findIC(x,Cin) 
  if x==0
    C0 = Cin;
  else
    C0 = 0.0;
  end
end 

[YeEn Kim]

Hi! Thank you again! But, I still get all the z values = 0 in the graph. I tried to debug it but seems like there's no problem...

[Torsten]

What are the calculated constants, especially Da, u0 and constant.

[YeEn Kim]

constant=1.0205333333333

Da=1.751055844870244e-10

u0=0.00194 (velocity not related to u in the function)

[Torsten]

If Da is that small, you'll have to set the final time of integration to a very high value to see anything developing. Set tend to 1e6, e.g.

Answer 2

[Bill Greene]

If you are solving a PDE with either cylindrical symmetry (m=1, your case) or spherical

symmetry (m=2), and your left boundary is at r=0, the symmetry condition requires that

the derivative of the dependent variable(s) (c in your case) equal zero. So pdepe simply ignores

the left boundary condition you specify in the boundary condition function and substitutes this

symmetry condition.

Since you have at the right end, no source term in your PDE, and initial conditions

equal zero, the solution is zero over the whole region for all time.

[Torsten]

@Bill Greene

Yes, that's true. I overlooked the m=1 setting.