Help with numerical solution to Heat & Mass diffusion PDEs

7 次查看（过去 30 天）
Julio
Julio 2022-1-26
评论： Torsten 2022-1-27
I'm trying to recreate a model from a paper, however I am unsure how to implement one of the differential terms. The paper is attached along with a more detailed description of my problem. I do not know how to implement the dX/dT term from the heat equation into the model,
%% Main
clear
close all
global counter
counter = 0;
x = 0.01;
deltax = 0.09*x ;
N = floor(x/deltax + 1 );
% y00= 700.6
y0(1:N,1) = 273+25 ; % initial Temperature
y0(N+1:N+N,1) = 0.4014 ; % initial moisture (wet basis (g water/g product))
y0(N+N+1:N*3,1) = 4.6315e-09/(12*12*2*10^(-6)*284*10^(-3)); % initial vapour g vapour/g product
tSpan = [0;80*60] ;
[tSol,ysol] = ode15s(@(t,y) NumericalCalc(t,y,x,y0), tSpan, y0);
TempSol = ysol(:,1:(N)) ;
MoistSol = ysol(:,N+1:N+N) ;
VapourSol = ysol(:,N+N+1:3*N);
plot(tSol,TempSol(:,1)-273)
figure()
plot(tSol,MoistSol(:,end))
figure()
plot(tSol,MoistSol(:,1))
%% NumericalCalc.m
function fval = NumericalCalc(t,y,x,y0)
%%% When referecing temperature at point "i", reference point is "i"
%%% When Referencing moisture at point "i", referecne point is "N+i"
%%% When Referencing vapour at point "i", reference point is N+N+i
global dXdtSave
global counter
counter = counter + 1;
N = length(y)/3 ; % number of nodes (dividing by 3 due to having moisture, vapour and temperature in y)
deltax = x*0.09;
%% Bread Parameters
% define dX/dt
dXdt= zeros(N*3,1);
if counter == 1
dXdtSave = dXdt;
end
Tinf = 273+210; % ambient temp
k = 0.07 ;
cp = 3500;
hc = 0.5;
Dw = 1.35*10^(-10);
Vair = 0.0205; %gram water/gram air
Temp = [0.01 2 4 10 14 18 20 25 30 34 40 44 50 54 60 70 80 90 96 100 110 120 130 140 150 160 180 200 220 240 260 280 300 320 340 360 370];
Psat = [ 0.006 0.007 0.008 0.0121 0.0158 0.0204 0.0231 0.0313 0.0419 0.0526 0.0729 0.0899 0.122 0.148 0.197 0.308 0.468 0.693 0.866 1.001 1.42 1.96 2.67 3.57 4.7 6.1 9.9 15.35 22.89 33.03 46.31 63.33 84.76 111.4 144.1 184.2 207.7];
SA = 0.12*0.02*2;
volume = 0.12*0.12*0.2;
a = 0.12/0.056;
b = a ;
a1 = 1+a^2;
b1=1+b^2;
ep = 0.9 ;
er = 0.95;
Fsp = 2/(pi*a*b)*(log(sqrt((a1*b1)/(1+a^2+b^2)+a*sqrt(b1)*atan(a/sqrt(b1))+b*sqrt(a1)*atan(b/sqrt(a1))-a*atan(a)-b*atan(b)))) ;
stf = 5.67*10^(-8);
VapourVolume = 0.12*0.12*deltax*0.7;
BreadVolume = 0.12*0.12*deltax ;
%% start iteration
for i=1:N
% update parameters
rho = 170+284*y(N+i) ;
gProduct = rho*BreadVolume*1000;
Dv = 9.0*10^(-12)*y(i)^2;
hv = (3.2*10^9)/(y(i)^3) ;
hw = 1.4*10^(-3)*y(i)+0.27*y(N+i)-4.0*10^(-4)*y(i)*-y(N+i)-0.77*y(N+i)^2 ;
r = -2.61E+03*(y(i)) + 2.51E+06; %J/kg
if i == 1
dXdt(N+i) = Dw*((y(N+i+2)-2*y(N+i+1)+y(N+i))/(deltax^2)); %liq water diffusion
dXdt(i) = k/(rho*cp)*((y(i+2)-2*y(i)+y(i))/(deltax^2))+1/cp*r*dXdtSave(N+i); %temp diffusion
%calculate water vapour, assuming saturation
[m,j] = min(abs(Temp-(y(i)-273)));
Pi = Psat(j);
n = Pi*VapourVolume/(8.205*10^(-5)*y(i)); %moles water in pore space
gramV = n*18;
y(N+N+i) = gramV/gProduct; %update vapour value
y(N+i) = y(N+i) - gramV/gProduct;%Subtract water vapour from liquid water
dXdt(N+N+i) = Dv*((y(N+N+i+2)-2*y(N+N+i+1)+y(N+N+i))/(deltax^2)); %Vap diffusion
%recalculate vapour as it may have changed. Change to Pressure to
%compare to saturation
gramV = y(N+N+i)*gProduct;
n = gramV/18;
Pi = n*8.205*10^(-5)*y(i)/(0.12*0.12*deltax*0.7);
[m,j] = min(abs(Temp-(y(i)-273)));
%determine if any vapour condenses
if Pi - Psat(j) <= 0
else
condenseP = Pi-Psat(j);
n = condenseP*VapourVolume/(8.205*10^(-5)*y(i));%mols of condensed
cWater=n*18; %grams of condensed water
y(N+i)=y(N+i)+ cWater/gProduct; %add condense water to liq water (units Kg Water/Kg Product or g/g)
y(N+N+i) = y(N+N+i)-cWater/gProduct; %subtract condensed water from vapour
end
% dXdt(N+i) = Dw*((y(N+i+2)-2*y(N+i+1)+y(N+i))/(deltax^2)); %liq water
elseif i== N
%Water
c = hw*(-y(N+i)+0); %%%NOTE REVERSED SIGNS as was gaining water
Wi = c*deltax+y(N+i-1);
dXdt(N+i) = Dw*((Wi-2*y(N+i-1)+y(N+i-2))/(deltax^2)); %liq water
hr = (stf*(Tinf^2+y(N)^2)*(Tinf+y(N)))/(1/ep+1/er-2+1/Fsp) ;
b = hc*(Tinf-y(N))+hr*(Tinf-y(N));
a = r*997*Dw*hw*(y(N+1+i)-0); % using B.C for dX/dr
Ti = ((b-a)/k)*deltax + y(i-1); %%% NOTE REVERSED SIGNS as was loosing heat (-k --> k)
dXdt(i) = k/(rho*cp)*(Ti-2*y(i-1)+y(i-2))/(deltax^2) + 1/cp*r*dXdtSave(N+i); %temp
%calculate water vapour, assuming saturation
[m,j] = min(abs(Temp-(y(i)-273)));
Pi = Psat(j);
n = Pi*VapourVolume/(8.205*10^(-5)*y(i)); %moles water in pore space
gramV = n*18;
y(N+N+i) = gramV/gProduct;%update vapour value
y(N+i) = y(N+i) - gramV/gProduct;%Subtract water vapour from liquid water
g = hv*(y(N+N+i)-Vair); % Reversed signs
Vi = g*deltax+y(N+N+i);
dXdt(N+N+i) = Dv*((Vi-2*y(N+N+i-1)+y(N+N+i-2))/(deltax^2)); %Vap
%Condense
gramV = y(N+N+i)*gProduct;
n = gramV/18;
Pi = n*8.205*10^(-5)*y(i)/VapourVolume;
[m,j] = min(abs(Temp-(y(i)-273)));
if Pi - Psat(j) <= 0
else
condenseP = Pi-Psat(j);
n = condenseP*VapourVolume/(8.205*10^(-5)*y(i));%mols of condensed
cWater=n*18; %grams of condensed water
y(N+i)=y(N+i)+ cWater/gProduct; %add condense water to liq water
y(N+N+i) = y(N+N+i)-cWater/gProduct; %subtract condensed water from vapour
end
else
dXdt(N+i) = Dw*(y(N+i+1)-2*y(N+i)+y(N+i-1))/(deltax^2); % liq water
dXdt(i) = k/(rho*cp)*(y(i+1)-2*y(i)+y(i-1))/(deltax^2)+1/cp*r*dXdtSave(N+i); %temp
%calculate water vapour
[m,j] = min(abs(Temp-(y(i)-273)));
Pi = Psat(j);
n = Pi*VapourVolume/(8.205*10^(-5)*y(i)); %moles water in pore space
gramV = n*18;
y(N+N+i) = gramV/gProduct; %update vapour
y(N+i) = y(N+i) - gramV/gProduct; %Subtract water vapour from liquid water
dXdt(N+N+i) = Dv*(y(N+N+i+1)-2*y(N+N+i)+y(N+N+i-1))/(deltax^2); % Vapour
%calculate condensed water
gramV = y(N+N+i)*gProduct;
n = gramV/18;
Pi = n*8.205*10^(-5)*y(i)/VapourVolume;
[m,j] = min(abs(Temp-(y(i)-273)));
if Pi - Psat(j) <= 0
else
condenseP = Pi-Psat(j);
n = condenseP*VapourVolume/(8.205*10^(-5)*y(i));%mols of condensed
cWater=n*18; %grams of condensed water
y(N+i)=y(N+i)+ cWater/gProduct; %add condense water to liq water
y(N+N+i) = y(N+N+i)-cWater/gProduct; %subtract condensed water from vapour
end
% dXdt(N+i) = Dw*(y(N+i+1)-2*y(N+i)+y(N+i-1))/(deltax^2); % liq water
end
dXdtSave = dXdt(:);
end
%Extract Fval from dXdt
fval = dXdt(:);
  3 个评论
显示 1更早的评论
Julio
Julio 2022-1-27
I am trying to avoid using pdepe, later on I will need to include shrinkage into the equations for my own work, I've tried including a moving mesh/changing the equations to include dimensionless length into pdepe before but had difficulties
While V doesn't appear in the other equations explicitly, it relies on the solution of the heat equation to calculate the vapour at each step.
Subbing looks like it works, thanks for that.
Other issues I'm having are that
  • I’ve had to change the signs (+’ve, -‘ve) in the boundary conditions in order to get the variable going in the correct direction.
  • The temperature does not plateau like it does in the paper, I think that the enthalpy of condensation needs to be considered between steps 4 and 5 although I don’t know how to implement that. (, I need an appropriate term with time)
  • The vapour content in the centre of the particle continues to increase past reasonable values,
  • I think somehow I need to test if there is any water available to be turned into vapour, I tried doing an if statement however it caused more problems.
Or should I move these questions to another forum as they are more at the understanding/conceptualization rather than coding level?
Torsten
Torsten 2022-1-27
At least to validate the results from your code, you should implement the above equations in pdepe for comparison.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Fluid Mechanics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by