Solving nonlinear advection diffusion equation with pdepe

I am trying to solve the following nonlinear advection diffusion equation with pdepe :

function [x,h] = Transient();
%% Initialization
theta = 50; % degrees - angle of repose of material
n = 8; % r.p.m. - rotation speed of tube
rho = 660; % kg/m3 - bulk density of material
R = 0.075; % m - inner radius of tube
beta = 3; % degrees - inclination angle of kiln
L_tube = 2; % m,lenght of the kiln
mass_flow_rate = 0.5; % kg/min - mass flow of material
x_points=100;
t_points=100;

%% Solving initil condition
xspan_1=[0 L_tube];
h0 = 0.000001; %  initial condition (bed height at the discharge end - average of particle diameter)
sol = ode45( @saeman_bed ,xspan_1 ,h0); %  solver for unsteady nonliner ode
x = linspace(0,L_tube,x_points);
u_0 = deval(sol,x);

%% Solving with pdepe
m = 0;
x = linspace(0,1,x_points); 
t = linspace(0,2,t_points); 
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
% Extract the first solution component as u.
u = sol(:,:,1)

function [c,f,s] = pdex1pde(x,t,u,DuDx)
f_H=power((2*u/R)-power(u,2)/power(R,2),0.5);
u_T=2*pi*n*R;

c = f_H;
f = (u_T*R*cotd(theta))*power(f_H,1.5).*DuDx;
s = u_T*tand(beta)/sind(theta)*f_H.*power(1-f_H,0.5).*DuDx;
end

% --------------------------------------------------------------

function u0 = pdex1ic(x)
u0 = transpose(u_0);
% -------------------------------
end

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul+0.000001 % exit 
ql = 0;
pr =(3*tan(deg2rad(theta))* mass_flow_rate/rho) / (4*pi*n*power(R,3)) * power((2*ul/R) - power((ul/R),2),(-3/2)) -tan(deg2rad(beta))/cos(deg2rad(theta));
qr = -1;
end

%% function for initial condition
function [dhdx,h] = saeman_bed (z,h);
dhdx = (3*tan(deg2rad(theta))* mass_flow_rate/rho) / (4*pi*n*power(R,3)) * power((2*h/R) - power((h/R),2),(-3/2)) -tan(deg2rad(beta))/cos(deg2rad(theta));
end
end

I have problems with boundary conditions "Unexpected output of BCFUN. For this problem BCFUN must return four column vectors of length 100." . Maybe somebody knows how properly set boundary conditions and pdex1pde function.

Thank you in advance!

Answer 1

[Torsten]

function [x,h] = Transient();

%% Initialization

theta = 50; % degrees - angle of repose of material
n = 8; % r.p.m. - rotation speed of tube
rho = 660; % kg/m3 - bulk density of material
R = 0.075; % m - inner radius of tube
beta = 3; % degrees - inclination angle of kiln
L_tube = 2; % m,lenght of the kiln
mass_flow_rate = 0.5; % kg/min - mass flow of material

x_points=100;
t_points=100;

%% Solving initial condition

xspan_1=[0 L_tube];
h0 = 0.000001; %  initial condition (bed height at the discharge end - average of particle diameter)
sol = ode45( @saeman_bed ,xspan_1 ,h0); %  solver for unsteady nonliner ode
x = linspace(0,L_tube,x_points);
u0 = deval(sol,x)

%% Solving with pdepe

m = 0;
x = linspace(0,1,x_points); 
t = linspace(0,2,t_points); 
sol = pdepe(m,@pdex1pde,@(x0)pdex1ic(x0,x,u0),@pdex1bc,x,t);

% Extract the first solution component as u.

u = sol(:,:,1)


function [c,f,s] = pdex1pde(x,t,u,DuDx)
f_H=power((2*u/R)-power(u,2)/power(R,2),0.5);
u_T=2*pi*n*R;

c = f_H;
f = (u_T*R*cotd(theta))*power(f_H,1.5).*DuDx;
s = u_T*tand(beta)/sind(theta)*f_H.*power(1-f_H,0.5).*DuDx;

end

% --------------------------------------------------------------

function u0 = pdex1ic(x0,x,u)
u0 = interp1(x,u,x0);
% -------------------------------
end


function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul+0.000001 % exit 
ql = 0;
pr =(3*tan(deg2rad(theta))* mass_flow_rate/rho) / (4*pi*n*power(R,3)) * power((2*ul/R) - power((ul/R),2),(-3/2)) -tan(deg2rad(beta))/cos(deg2rad(theta));
qr = -1;
end


%% function for initial condition

function [dhdx,h] = saeman_bed (z,h);
dhdx = (3*tan(deg2rad(theta))* mass_flow_rate/rho) / (4*pi*n*power(R,3)) * power((2*h/R) - power((h/R),2),(-3/2)) -tan(deg2rad(beta))/cos(deg2rad(theta));
end

end

[Alina Mexdi]

Dear Torsten,

Thank you for your help butnow I have a warning " Warning: Failure at t=2.870379e-04. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (8.673617e-19) at time t."

[Torsten]

Your boundary condition at the right end is incorrect.

Note that pr and qr are to be set such that

pr + qr*f = 0

not

pr + qr*dh/dz=0

[Alina Mexdi]

I am new in matlab pdepe function so seems I understand in wrongly, I am trying to code right boundary with this expression :

Can you explain me please how it could be done in a proper way ?

Do I underatsnd correctly that instead of :

qr=-1

it should be

qr = -(u_T*R*cotd(theta))*power(f_H,1.5)

[Torsten]

No,

qr = -1/(u_T*R*cotd(theta))*power(f_H,1.5))

And you use "ul" in the calculation of "pr". It has to be "ur".

And where is the 1/3 in the calculation of f ?

And you set h(t,z=0) = -0.000001 as boundary condition. Is this negative value for h correct ?

[Alina Mexdi]

Thank you very much Mr. Torsten . Now it converges but if I plot the results is shows me the following:

Meaning that nothing changes with time. But in my case it should show me evolution of filling degre of bed in the rotary kiln.

I have tried another BC:

pl = ul-0.00008;
ql = ones(x_points,1).*0;  
pr =(3*tan(deg2rad(theta))* mass_flow_rate/rho) / (4*pi*n*power(R,3)) * power((2*ur/R) - power((ur/R),2),(-3/2)) -tan(deg2rad(beta))/cos(deg2rad(theta));
qr = ones(x_points,1).*-1;

ql and qr are set up in this way because otherwise it says that "Unexpected output of BCFUN. For this problem BCFUN must return four column vectors of length 50."

and the resulting plot i:

which is exactly what i would like to see as a result. I am doing something crazy or totally wrong? Or why can't I get right evoulution of the function with respect to time with right BC?

[Torsten]

Please include again the complete code you are using at the moment.

[Alina Mexdi]

This is the code which produces kind of "right" evolution with time but as you noted wrong BC :

function [x,h] = Transient();
%% Initialization
theta = 36.44; % degrees - angle of repose of material
n = 8; % r.p.m. - rotation speed of tube
rho = 660; % kg/m3 - bulk density of material
R = 0.075; % m - inner radius of tube
beta = 3; % degrees - inclination angle of kiln
L_tube = 2; % m,lenght of the kiln
Time=20;
mass_flow_rate = 0.5; % kg/min - mass flow of material
x_points=50;
t_points=100;

%% Solving initil condition
xspan_1=[0 L_tube];
h0 = 0.00008; %  initial condition (bed height at the discharge end - average of particle diameter)
sol = ode45( @saeman_bed ,xspan_1 ,h0); %  solver for unsteady nonliner ode
x = linspace(0,L_tube,x_points);
u_0 = deval(sol,x)

%% Solving with pdepe

m = 0;
x = linspace(0,L_tube,x_points); 
t = linspace(0,Time,t_points); 
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);

% Extract the first solution component as u.
u = sol(:,:,1)

% A surface plot is often a good way to study a solution.
surf(x,t,u) 
title('Numerical solution.')
xlabel('Distance x')
ylabel('Time t')

% A solution profile can also be illuminating.
figure
plot(x,u(end,:))
title('Solution at t = final')
xlabel('Distance x')
ylabel('Bed height')

function [c,f,s] = pdex1pde(x,t,u,DuDx)
f_H=power((2*u/R)-power(u,2)/power(R,2),0.5);
u_T=2*pi*n*R;

c = f_H;
f = ((u_T/3)*R*cotd(theta))*power(f_H,1.5).*DuDx;
s = u_T*tand(beta)/sind(theta)*f_H.*power(1-f_H,0.5).*DuDx;

end
% --------------------------------------------------------------
function u0 = pdex1ic(x);
u0 = transpose(u_0);
% -------------------------------
end

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t);
    
pl = ul-0.00008;
ql = ones(x_points,1).*0;
pr =(3*tan(deg2rad(theta))* mass_flow_rate/rho) / (4*pi*n*power(R,3)) * power((2*ur/R) - power((ur/R),2),(-3/2)) -tan(deg2rad(beta))/cos(deg2rad(theta));
qr = ones(x_points,1).*-1;
end


%% function for initial condition
function [dhdx,h] = saeman_bed (z,h);
dhdx = (3*tan(deg2rad(theta))* mass_flow_rate/rho) / (4*pi*n*power(R,3)) * power((2*h/R) - power((h/R),2),(-3/2)) -tan(deg2rad(beta))/cos(deg2rad(theta));
end
end

And this is the code with your corrections and right BC:

function [x,h] = Transient();
%% Initialization
theta = 36.44; % degrees - angle of repose of material
n = 8; % r.p.m. - rotation speed of tube
rho = 660; % kg/m3 - bulk density of material
R = 0.075; % m - inner radius of tube
beta = 3; % degrees - inclination angle of kiln
L_tube = 2; % m,lenght of the kiln
mass_flow_rate = 0.5; % kg/min - mass flow of material
x_points=50;
t_points=30;

%% Solving initial condition
xspan_1=[0 L_tube];
h0 = 0.000001; %  initial condition (bed height at the discharge end - average of particle diameter)
sol = ode45( @saeman_bed ,xspan_1 ,h0); %  solver for unsteady nonliner ode
x = linspace(0,L_tube,x_points);
u0 = deval(sol,x)

%% Solving with pdepe
m = 0;
x = linspace(0,2,x_points); 
t = linspace(0,5,t_points); 
sol = pdepe(m,@pdex1pde,@(x0)pdex1ic(x0,x,u0),@pdex1bc,x,t);
% Extract the first solution component as u.
u = sol(:,:,1)
% A surface plot is often a good way to study a solution.
surf(x,t,u) 
title('Numerical solution')
xlabel('Distance x')
ylabel('Time t')

% A solution profile can also be illuminating.
figure
plot(x,u(end,:))
title('Solution at t = final')
xlabel('Distance x')
ylabel('Bed height')

function [c,f,s] = pdex1pde(x,t,u,DuDx)
f_H=power((2*u/R)-power(u,2)/power(R,2),0.5);
u_T=2*pi*n*R;
c = f_H;
f = ((u_T/3)*R*cotd(theta))*power(f_H,1.5).*DuDx;
s = u_T*tand(beta)/sind(theta)*f_H.*power(1-f_H,0.5).*DuDx;
end
% --------------------------------------------------------------
function u0 = pdex1ic(x0,x,u)
u0 = interp1(x,u,x0);

% -------------------------------
end
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul-0.000001; % exit 
ql = 0;
pr =(3*tan(deg2rad(theta))* mass_flow_rate/rho) / (4*pi*n*power(R,3)) * power((2*ur/R) - power((ur/R),2),(-3/2)) -tan(deg2rad(beta))/cos(deg2rad(theta));
qr = -1/(2*pi*n*R*R*cotd(theta))*power(power((2*ur/R)-power(ur,2)/power(R,2),0.5),1.5);
end
%% function for initial condition
function [dhdx,h] = saeman_bed (z,h);
dhdx = (3*tan(deg2rad(theta))* mass_flow_rate/rho) / (4*pi*n*power(R,3)) * power((2*h/R) - power((h/R),2),(-3/2)) -tan(deg2rad(beta))/cos(deg2rad(theta));
end
end

P.S. in this case I don't understand what is passed to pdex1ic(x0,x,u0) ,namely what is x0 ,becauseit isnot defined anyway.

[Torsten]

%% Initialization

theta = 36.44; % degrees - angle of repose of material
n = 8; % r.p.m. - rotation speed of tube
rho = 660; % kg/m3 - bulk density of material
R = 0.075; % m - inner radius of tube
beta = 3; % degrees - inclination angle of kiln
L_tube = 2; % m,lenght of the kiln
mass_flow_rate = 0.5; % kg/min - mass flow of material
x_points=50;
t_points=30;

%% Solving initial condition

xspan_1=[0 L_tube];
h0 = 0.000001; %  initial condition (bed height at the discharge end - average of particle diameter)
sol = ode45( @saeman_bed ,xspan_1 ,h0); %  solver for unsteady nonliner ode
x = linspace(0,L_tube,x_points);
u0 = deval(sol,x);


%% Solving with pdepe

m = 0;
x = linspace(0,2,x_points); 
t = linspace(0,5,t_points); 
sol = pdepe(m,@pdex1pde,@(x0)pdex1ic(x0,x,u0),@pdex1bc,x,t);

% Extract the first solution component as u.

u = sol(:,:,1)

% A surface plot is often a good way to study a solution.

surf(x,t,u) 
title('Numerical solution')
xlabel('Distance x')
ylabel('Time t')



% A solution profile can also be illuminating.

figure
plot(x,u(1,:),x,u(end,:))
title('Solution at t = final')
xlabel('Distance x')
ylabel('Bed height')



function [c,f,s] = pdex1pde(x,t,u,DuDx)

f_H=power((2*u/R)-power(u,2)/power(R,2),0.5);
u_T=2*pi*n*R;

c = f_H;
f = ((u_T/3)*R*cotd(theta))*power(f_H,1.5).*DuDx;
s = u_T*tand(beta)/sind(theta)*f_H.*power(1-f_H,0.5).*DuDx;

end

% --------------------------------------------------------------

function u0 = pdex1ic(x0,x,u)

u0 = interp1(x,u,x0);

% -------------------------------

end

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)

pl = ul-0.000001; % exit 
ql = 0;
pr =(3*tan(deg2rad(theta))* mass_flow_rate/rho) / (4*pi*n*power(R,3)) * power((2*ur/R) - power((ur/R),2),(-3/2)) -tan(deg2rad(beta))/cos(deg2rad(theta));
qr = -1/(2*pi*n*R/3*R*cotd(theta))*power(power((2*ur/R)-power(ur,2)/power(R,2),0.5),1.5);

end

%% function for initial condition

function [dhdx,h] = saeman_bed (z,h);

dhdx = (3*tan(deg2rad(theta))* mass_flow_rate/rho) / (4*pi*n*power(R,3)) * power((2*h/R) - power((h/R),2),(-3/2)) -tan(deg2rad(beta))/cos(deg2rad(theta));

end

end

[Alina Mexdi]

Dear Torsten,

Still as can be seen the curve doesn't change with the time

Probably I am doing something wrong because my initial profile u (x,0) actually should be considered as a final steady state solution and by transient regime I want to see the evolution of this fucntion till that steady state moment. like here:

[Torsten]

But if you start with a steady-state profile obtained from ODE45, why do you expect it should change in the transient calculation ?

In the plot where you see an evolution of the profile, you start with a profile that is not the one obtained from ODE45 as you can see when you compare the two surfaces at time t=0. So you will see an evolution as it is not steady-state.

[Alina Mexdi]

Thank you very much , now verything became clear. The last question please , I am calculating steady state solution with ode45, and I'm calculating transient one with pdepe and pdepe slightly underestimate the function :

Is it a feature/drawback/difference of the solvers , or I should look carefully for mistakesin my code ?

[Torsten]

I don't understand why the ODE45 solution is the steady-state solution of the problem.

The steady-state solution is the solution of the transient problem if you neglect time-dependent terms. In the present problem, remove F_h^1/2 * dh/dt from your PDE equation and leave the boundary conditions as they are. You'll arrive at a second-order ODE for h with one boundary condition at the left and one boundary condition at the right endpoint of the spatial interval. This could be solved using BVP4C.

Maybe the problem you solve with ODE45 is an approximation to the steady-state profile.