function [x,h] = Transient()
theta = 50;
n = 8;
rho = 660;
R = 0.075;
beta = 3;
L_tube = 2;
mass_flow_rate = 0.5;
x_points=100;
t_points=100;
xspan_1=[0 L_tube];
h0 = 0.000001;
sol = ode45( @saeman_bed ,xspan_1 ,h0);
x = linspace(0,L_tube,x_points);
u0 = deval(sol,x)
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);
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
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 [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
8 Comments
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/437971-the-following-code-is-not-working#comment_656275
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/437971-the-following-code-is-not-working#comment_656275
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/437971-the-following-code-is-not-working#comment_656284
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/437971-the-following-code-is-not-working#comment_656284
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/437971-the-following-code-is-not-working#comment_656322
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/437971-the-following-code-is-not-working#comment_656322
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/437971-the-following-code-is-not-working#comment_656410
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/437971-the-following-code-is-not-working#comment_656410
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/437971-the-following-code-is-not-working#comment_656461
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/437971-the-following-code-is-not-working#comment_656461
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/437971-the-following-code-is-not-working#comment_656617
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/437971-the-following-code-is-not-working#comment_656617
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/437971-the-following-code-is-not-working#comment_656647
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/437971-the-following-code-is-not-working#comment_656647
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/437971-the-following-code-is-not-working#comment_656757
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/437971-the-following-code-is-not-working#comment_656757
Sign in to comment.