pdepe help! The solution gives 0...

3 views (last 30 days)
YeEn Kim
YeEn Kim on 27 Mar 2022
Commented: Torsten on 2 Apr 2022
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

Accepted Answer

Torsten
Torsten on 27 Mar 2022
IC = @(x) findIC(x,Cin);
...
function C0=findIC(x,Cin)
if x==0
C0 = Cin;
else
C0 = 0.0;
end
end
  4 Comments
YeEn Kim
YeEn Kim on 27 Mar 2022
constant=1.0205333333333
Da=1.751055844870244e-10
u0=0.00194 (velocity not related to u in the function)
Torsten
Torsten on 27 Mar 2022
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.

Sign in to comment.

More Answers (1)

Bill Greene
Bill Greene on 2 Apr 2022
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.

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!