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
采纳的回答
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
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Antennas and Electromagnetic Propagation 的更多信息
Translated by 
