PDEPE coupled boundary condition help
1 view (last 30 days)
Show older comments
Hi there,
I am trying to simulate heat and mass transfer during the drying of a biomass particle (spherical), such that I can have plots of Temperature and Moisture content. I have used pdepe to try and achieve this. I have had it working with simplier cases but currently I am having difficulty with one of my boundary conditions. Code lines which have *** next to them are my attempted solution. I know global variables are bad but I was just trying to see if this would work.
I have attached a file with the equations I am trying to solve with the initial and boundary conditions.
%Main -------------------
clear
close all
global prev_ur; %***
global prev_t; %***
prev_ur = [293;3.5] %***
prev_t = 0.0001; %***
xmesh = linspace(0,0.005,50);
tspan = linspace(0,300,50);
m = 2;
sol = pdepe (m,@pdefun,@pdeic,@pdebc,xmesh,tspan);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
%pdefun---------------------
function [c,f,s] = pdefun(x,t,u,dudx)
rho = 11.35*u(2)^4-92.25*u(2)^3+277.04*u(2)^2-402.65*u(2)+1110.60 ;
cp = -44.93*u(2)^4+382.13*u(2)^3-1246.3*u(2)^2+2123.60*u(2)+1678.10;
ktherm = -4.70*10^(-3)*u(2)^4+3.94*10^(-3)*u(2)^3-1.25*10^(-2)*u(2)^2+0.21*u(2)+0.16 ;
c = [(rho*cp)/ktherm; 1/(7.4*10^(-9))];
f = [1; 1] .* dudx;
s = [0;0] ;
end
%pdeic---------------------
function u0 = pdeic(x)
u0 = [(273+20); 3.5];
end
%pdebc-----------------------------
function [pl, ql, pr, qr] = pdebc(xl,ul,xr,ur,t)
den = 11.35*ur(2)^4-92.25*ur(2)^3+277.04*ur(2)^2-402.65*ur(2)+1110.60 ;
cp = -44.93*ur(2)^4+382.13*ur(2)^3-1246.3*ur(2)^2+2123.60*ur(2)+1678.10;
ktherm = -4.70*10^(-3)*ur(2)^4+3.94*10^(-3)*ur(2)^3-1.25*10^(-2)*ur(2)^2+0.21*ur(2)+0.16 ;
global prev_ur %***
global prev_t %***
diffx = ur(2)-prev_ur(2); %***
deltar = 0.005;
difft = t-prev_t;%***
prev_ur = ur
pl = [0; 0];
ql = [1; 1];
pr = [-(11.72*(ur(1)-443))/(ktherm)-(den*deltar*(diffx/difft)*(2256400+1880*(ur(1)-443)))/(ktherm); (1.63*10^(-5))/(7.4*10^(-9))*(ur(2)-0.0032)];
qr = [-1;1];
%Note: code works if using pr = [-(11.72*(ur(1)-443))/(ktherm); (1.63*10^(-5))/(7.4*10^(-9))*(ur(2)-0.0032)]; however not accurate compared to experimental
end
Thank you for your help
0 Comments
Accepted Answer
Torsten
on 24 May 2021
Edited: Torsten
on 24 May 2021
dX^bar/dt = 3/R * D_ap * dX/dr (@r=R)
For dX/dr (@r=R) you can insert the right-hand side of your second boundary condition for X.
The boundary condition for T can't depend on a numerical value, namely deltar. Or what does this deltar stand for else than R/(number of discretization points in r-direction -1) ? Could you elaborate ?
Further, the change of the radius of the particle should be accounted for, I guess.
7 Comments
Torsten
on 27 May 2021
The coordinate for r at each time t runs between r = 0 and r = R(t).
If you normalize this r-coordinate as r' = r/R(t), the r' coordinate runs between r' = 0/R(t) = 0 and r' = R(t)/R(t) = 1.
More Answers (0)
See Also
Categories
Find more on Eigenvalue Problems in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!