# PDEPE coupled boundary condition help

9 views (last 30 days)
julian nylen on 24 May 2021
Commented: Torsten on 27 May 2021
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

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.
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.