PDEPE coupled boundary condition help

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

Answer 1

[Torsten]

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.

[Julio]

Thank you for your quick reply Torsten. Would you mind telling me how you derived your dX^bar/dt equation? The model is looking better, although still not great, most likely due to that deltar. I hadn't thought about the equations that much, but thank you for pointing out that B.C for T can't depend on a numerical value.

Admittally these equations are not mine, I am trying to model the data found in this paper https://doi.org/10.1016/j.applthermaleng.2012.04.025.

The authors unfortunetly don't explicitly say what deltar is, although your assumption sounds good. I should have elborated that X^bar is the medium moisture content along the section. My assumption was that deltar would then be the "r" for that section/ring. If you could make a guess as to what they mean from the paper, I'd be very thankful.

Again you are correct that the change in particle size due to shrinkage should be considered. However the authors of the above paper did not consider it. You might find another question from me regarding that phenomena when I begin doing models for my biomass :).

Another question when specifying the length in xmesh, is this length the diameter or radius of the object?

Again thank you for your help.

[Julio]

I just modelled the data again ignoring the dX^bar/dt side of the equation and found that it does fit the experimental data quite well (I must have made a mistake somewhere when I modelled it the first time). Whether my assumption or yours is correct about it I don't think it matters. If we assume that we have a large number of nodal points/sections/rings then deltar approaches zero, effectively ignoring that half of the RHS.

Regardless, I would still be interested in how you derived the dX^bar/dt equation. (I also realised that the length in xmesh has to be radius for my case)

Again thank you for your help

[Torsten]

The rate of change of mean moisture content within the particle is equal to the flux across its boundaries, divided by the volume of the particle.

Flux across the boundaries:

r=0 : flux_0 = 0

r=R: flux_R = 4*pi*R^2 * D_ap * dX/dr (@r=R)

Thus dX^bar/dt = (flux_0 + flux_R)/(4/3*pi*R^3).

The values you specify in xmesh are the radii, thus in your code from above, the radius of the particle is assumed to be 0.005 [m].

[Julio]

Thank you for the derivation Torsten.

Regarding the change in radius due to shrinkage. From the litature it seems most common for the change to be represent by:

R = R_0 * aX/X_0

With R being the new radius, R_0 being the initial raduis, "a" being some correlation coefficeint and X average moisture content across the particle.

Is it possible to implmement this condition and solve with pdepe?

[Torsten]

The mesh for pdepe must remain fixed during the calculation.

So the only way is to introduce a dimensionless radius

r^~ = r/R(t)

and to make the corresponding coordinate transformation in your equations.

This way, the mesh always is on [0,1].

To get the average moisture content X^bar and dR/dt (which are needed in the transformed PDE) will also be a problem because in pdefun, you only have the value for X in one grid point, not over the complete radius. So, if possible at all, you will need quite a few tricks to make your modified model fit into the PDE scheme.

I think it will be much faster to discretize the equations of your new model and solve them using an ODE integrator with the method of lines.

[Julio]

Thanks for the suggestion Torsten,

Just one thing shouldn't r' = R(t)/r? In the shrinkage case R(t) will be smaller than r when t>0. Thus R(t)/r will give a value between 0 and 1

Or am I missinterpreting what you're saying.

[Torsten]

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.