pdepe function doesnt stop compiling

I'm trying to simulate the fluid drainage of a bubble column for three different bubble radii R with the help of the pdepe-function. It works if it only needs to simulate the process for one radius but as soon as I choose to do it with more it doesnt stop compiling.

Here is my code thus far:

par = getParameter();


m = 0;
xmesh = [0:par.dh: par.H];
tspan = [0: par.dt: par.Tend];

timeSteps = [0:0.05*par.Tend:par.Tend];


sol = pdepe(m, @model,@icfun, @bcfun, xmesh, tspan);

function par = getParameter()
%geometry
par.H  = 0.1;    % m
par.dh = 0.01;  % m 

% simulation time
par.Tend = 100;     %s
par.dt   = 0.01;  %s

%material

par.R1 = 0.0005;
par.R2 = 0.001;
par.R3 = 0.005;
par.rho = 1000; % kg/m^3
par.gamma = 0.04; % N/m
par.eta = 0.001; % Pa*s

%constants
par.k1 = 0.0066;
par.k2 = 0.161^(1/2);
par.g = 9.81; % m/s^2
end


function [c, f, s] = model(x,t,u,dudx)

par = getParameter();

c(1,1) = 1;
c(2,1) = 1;
c(3,1) = 1;

% flux term

f(1,1) = -((par.k1*par.R1^2)*par.eta^-1)*(par.rho*par.g*u(1)^2-...
    par.k2*(par.gamma*(2*par.R1)^-1)*u(1)^(1/2)*dudx(1));

f(2,1) = -((par.k1*par.R2^2)*par.eta^-1)*(par.rho*par.g*u(2)^2-...
    par.k2*(par.gamma*(2*par.R2)^-1)*u(2)^(1/2)*dudx(2));

f(3,1) = -((par.k1*par.R3^2)*par.eta^-1)*(par.rho*par.g*u(3)^2-...
    par.k2*(par.gamma*(2*par.R3)^-1)*u(3)^(1/2)*dudx(3));
s(1,1) = 0;
s(2,1) = 0;
 s(3,1) = 0;

end


% initial condition
function [u0]  = icfun(x)

par = getParameter();

u0(1,1) = 0.36;
u0(2,1) = 0.36;
u0(3,1) = 0.36;

end


% boundary conditions
function [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)
par = getParameter();

pl(1,1) = 0;
pl(2,1) = 0;
pl(3,1) = 0;
ql(1,1) = 1;
ql(2,1) = 1;
ql(3,1) = 1;
pr(1,1) = ur(1)-0.36;
pr(2,1) = ur(2)-0.36;
pr(3,1) = ur(3)-0.36;
qr(1,1) = 0;
qr(2,1) = 0;
qr(3,1) = 0;


end

Answer 1

[Torsten]

You set

par.rho*par.g*u(1)^2-par.k2*(par.gamma*(2*par.R1)^-1)*u(1)^(1/2)*dudx(1) = 0

at x=0. Are you sure this is correct ?

par = getParameter();


m = 0;
xmesh = [0:par.dh: par.H];
tspan = [0: par.dt: par.Tend];

timeSteps = [0:0.05*par.Tend:par.Tend];


sol = pdepe(m, @model,@icfun, @bcfun, xmesh, tspan);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3);
plot(xmesh,[u1(end,:);u2(end,:);u3(end,:)])

function par = getParameter()
%geometry
par.H  = 0.1;    % m
par.dh = 0.0001;  % m 

% simulation time
par.Tend = 100;     %s
par.dt   = 0.01;  %s

%material

par.R1 = 0.0005;
par.R2 = 0.001;
par.R3 = 0.005;
par.rho = 1000; % kg/m^3
par.gamma = 0.04; % N/m
par.eta = 0.001; % Pa*s

%constants
par.k1 = 0.0066;
par.k2 = 0.161^(1/2);
par.g = 9.81; % m/s^2
end


function [c, f, s] = model(x,t,u,dudx)

par = getParameter();

c(1,1) = 1;
c(2,1) = 1;
c(3,1) = 1;

% flux term

f(1,1) = -((par.k1*par.R1^2)*par.eta^-1)*(par.rho*par.g*u(1)^2-...
    par.k2*(par.gamma*(2*par.R1)^-1)*u(1)^(1/2)*dudx(1));

f(2,1) = -((par.k1*par.R2^2)*par.eta^-1)*(par.rho*par.g*u(2)^2-...
    par.k2*(par.gamma*(2*par.R2)^-1)*u(2)^(1/2)*dudx(2));

f(3,1) = -((par.k1*par.R3^2)*par.eta^-1)*(par.rho*par.g*u(3)^2-...
    par.k2*(par.gamma*(2*par.R3)^-1)*u(3)^(1/2)*dudx(3));
s(1,1) = 0;
s(2,1) = 0;
s(3,1) = 0;

end


% initial condition
function [u0]  = icfun(x)

par = getParameter();

u0(1,1) = 0.36;
u0(2,1) = 0.36;
u0(3,1) = 0.36;

end


% boundary conditions
function [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)
par = getParameter();

pl(1,1) = 0;
pl(2,1) = 0;
pl(3,1) = 0;
ql(1,1) = 1;
ql(2,1) = 1;
ql(3,1) = 1;
pr(1,1) = ur(1)-0.36;
pr(2,1) = ur(2)-0.36;
pr(3,1) = ur(3)-0.36;
qr(1,1) = 0;
qr(2,1) = 0;
qr(3,1) = 0;


end

[cloutdrop]

Not conciously at least. Could you explain what you mean please?

[Torsten]

This is your boundary condition at x=0:

pl(i,1) = 0  
ql(i,1) = 1 

means

pl(i,1) + ql(i,1) * f(i,1) = 0

, thus

0 + 1*-((par.k1*par.Ri^2)*par.eta^-1)*(par.rho*par.g*u(i)^2-par.k2*(par.gamma*(2*par.Ri)^-1)*u(i)^(1/2)*dudx(i)) = 0

[cloutdrop]

Then yes it is correct. It is stated that the liquid flux is zero at z = 0 which means that f=0 if i'm not mistaken.