PDEPE solver failure (spatial discretization) ; partial differential/derivative temperature modelling

Attempting to model temperature profile of an adsorption column based on concentration profile since this determines the heat of adsorption.

This is the function:

function[c,f,s] = pdefunEB(x,t,u,dudx)

global qmax_phenol K_L_phenol kf De u_s Rep Dp epsilon rho_liq cp_AC cp_feed H_ads h alpha

C = u(1) ;
q = u(2) ;
Tf = u(3) ;
dCdx = dudx(1) ;
dqdx = dudx(2) ;
dTdx = dudx(3) ;

Ts = 298.15 ; % adsorbent temp Kelvin

qstar = ((qmax_phenol)*(K_L_phenol)*C)/(1+(K_L_phenol)*C) ;
dqdt = kf*(qstar-q) ;

c = ones(3,1) ;
f = [De*dCdx;
        0;
        0];
s = [-u_s*dCdx-(1-epsilon)/epsilon*rho_liq*dqdt;
        dqdt;
        -u_s*dTdx+(((h*alpha*(Tf-Ts)-(H_ads/cp_AC))+(H_ads/cp_feed))*(dqdt*((1-epsilon)/epsilon)))] ;
end

%% initial condition
function u0 = icfunEB(x)
u0 = [0;
        0;
        298.15] ;
end 

%% boundary condition
function [p1,q1,pr,qr] = bcfunEB(x1,u1,xr,ur,t)

global C_phenol

C0 = C_phenol ;
p1 = [u1(1)-C0;0;-298.15] ;
q1 = [0;1;0] ;
pr = [0;0;0] ;
qr = [1;1;1] ;
end

this is where the error appears (Error using pdepe Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative)

mEB = 0 ;
sol1 = pdepe(mEB,@pdefunEB,@icfunEB,@bcfunEB,x,t) ; 
CEB = sol1(:,:,1) ; 
qEB = sol1(:,:,2) ;
TEB = sol1(:,:,3) ;

figure(2)
T0 = 298.15 ;
Cp1 = C1(:,end)./C01 ;
plot(t,TEB,'b-','LineWidth',2)
grid on
title('Column Temperature')
xlabel('time in minutes')
ylabel('temp in kelvin')

Heres a picture of my mathematical working before input into matlab:

Answer 1

[Torsten]

Please include the constants defined as globals so that we can test your code.

One obvious error is the boundary condition for T in the left boundary point. Use

p1 = [u1(1)-C0;0;u1(3)-298.15] ;

instead of

p1 = [u1(1)-C0;0;-298.15] ;

Another problem can be that f = 0 for equations (2) and (3). "pdepe" - as the name indicates - is designed for parabolic-elliptic partial differential equations (thus for equations with second-order spatial derivatives present). Your second equation is a simple ODE, your third equation is hyperbolic in nature.

[Josh]

global qmax_5HMF qmax_furfural qmax_acetic qmax_phenol 
global K_L_5HMF K_L_furfural K_L_acetic K_L_phenol
global C_5HMF C_furfural C_acetic C_phenol
global kf De u_s Rep Dp epsilon rho_liq cp_AC cp_feed H_ads h alpha

Dp = 2.5e-3 ; %% 4e-4 < Dp < 2.5e-3
epsilon = 0.427 ; %0.427
cp_AC = 0.751 ; %(Kutub Uddin source)
cp_feed = 4.177 ; %(see excel calcs)
H_ads = 33.15 ; %(Mingling Ren source)

F = 1.44120296 ; % inlet flow rate ; see excel 

u_s = 0.06 ; %% 0.06 < u_s < 0.24
A = F/u_s ;
L = 100*u_s ;
V = L*A ;
Dc = 2*(sqrt(A/pi)) ;

rho_AC = 1480 ; %% [kg/m3]
m_AC = epsilon*V*rho_AC ; %% 1480 kg/m3 AC

h = 0.017 ; %(Wang 2003)
specificA = 1000e3 ; %% 913e3 - 1413e3
alpha = specificA*m_AC ;

L/u_s %% L/u > 100

ans = 100

A*u_s %% A*u > F

ans = 1.4412

%% mass fractions
x_water = 0.894250861 ;
x_glucose = 0.025266561 ;
x_xylose = 0.037260453 ;
x_5HMF = 0.000371 ;
x_furfural = 0.00302 ;
x_acetic = 0.00634 ;
x_phenolic = 0.0335 ;

%% viscosities
mu_water = 0.00089002 ;
mu_glucose = 0.00124 ;
mu_xylose = 0.0000010 ;
mu_5HMF = 0.00092 ;
mu_furfural = 0.001494 ;
mu_acetic = 0.001037 ;
mu_phenolic = 0.0017844 ;

%% particle reynolds 
rho_liq = 1027.839426 ;
mu_liq = ((mu_water^x_water)+(mu_glucose^x_glucose)+(mu_xylose^x_xylose)+(mu_5HMF^x_5HMF)+(mu_furfural^x_furfural)+(mu_acetic^x_acetic)+(mu_phenolic^x_phenolic)) ;
Rep = (Dp*(u_s/60)*rho_liq)/(mu_liq) 

Rep = 4.9528e-04

%% pressure drop
deltaP = ((150/Rep)*(1-epsilon)+1.75)*(((1-epsilon)/(epsilon^3))*(L/Dp)*rho_liq*(u_s^2)) ;
deltaP %% deltaP < 20

deltaP = 1.1343e+10

kf = (0.0011*Rep)+0.0009 ; % film MTC [1/min]
De = ((1e-11)*Rep)+(1e-10) ; % diffusive term (Pergamon Source) [m2/min]

%% langmuir coefficients [m3/kg]
K_L_5HMF = 0.047 ; %(rajoriya source)
K_L_furfural = 0.047 ; %(rajoriya source)
K_L_acetic = 0.071889 ; %(Mutasim source)
K_L_phenol = 0.0154 ; %(xie source)

%% adsoprtion capacities [kg/kg]
qmax_5HMF = 86.21 ; %(rajoriya source)
qmax_furfural = 86.21 ; %(rajoriya source)
qmax_acetic =  108.5 ; %(wu source)
qmax_phenol = 214.13 ; %(xie source)

%% concentrations [kg/m3]
C_5HMF = 0.38122 ;
C_furfural = 3.10996 ;
C_acetic = 6.52088 ;
C_phenol = 34.413 ;

%%

x = linspace(0,L,201) ;
t = 0:10000:1000000 ; 

%% solve
m1 = 0 ;
sol1 = pdepe(m1,@pdefun1,@icfun1,@bcfun1,x,t) ; 
C1 = sol1(:,:,1) ; 
q1 = sol1(:,:,2) ;

m2 = 0 ;
sol2 = pdepe(m2,@pdefun2,@icfun2,@bcfun2,x,t) ; 
C2 = sol2(:,:,1) ; 
q2 = sol2(:,:,2) ;

m3 = 0 ;
sol3 = pdepe(m3,@pdefun3,@icfun3,@bcfun3,x,t) ; 
C3 = sol3(:,:,1) ; 
q3 = sol3(:,:,2) ;

m4 = 0 ;
sol4 = pdepe(m4,@pdefun4,@icfun4,@bcfun4,x,t) ; 
C4 = sol4(:,:,1) ; 
q4 = sol4(:,:,2) ;

%% plotting breakthrough curve
figure(1)
C01 = C_5HMF ;
Cp1 = C1(:,end)./C01 ;
plot(t,Cp1,'m-','LineWidth',1)
grid on
title('Breakthrough curve')
xlabel('time in minutes')
ylabel('C/C0')

hold on
C02 = C_furfural ;
Cp2 = C2(:,end)./C02 ;
plot(t,Cp2,'g--','LineWidth',1)

hold on
C03 = C_acetic ;
Cp3 = C3(:,end)./C03 ;
plot(t,Cp3,'b:','LineWidth',1)

hold on
C04 = C_phenol ;
Cp4 = C4(:,end)./C04 ;
plot(t,Cp4,'r-.','LineWidth',1)

%% ENERGY BAL

mEB = 0 ;
t=linspace(0,1.4,10);
sol1 = pdepe(mEB,@pdefunEB,@icfunEB,@bcfunEB,x,t) ; 
CEB = sol1(:,:,1) ; 
qEB = sol1(:,:,2) ;
TEB = sol1(:,:,3) ;

figure(2)
T0 = 298.15 ;
Cp1 = C1(:,end)./C01 ;
plot(t,TEB,'b-','LineWidth',2)
grid on
title('Column Temperature')
xlabel('time in minutes')
ylabel('temp in kelvin')

function[c,f,s] = pdefun1(x,t,u,dudx)

global qmax_5HMF K_L_5HMF kf De u_s Rep Dp epsilon rho_liq

C = u(1) ;
q = u(2) ;
dCdx = dudx(1) ;
dqdx = dudx(2) ;

qstar = ((qmax_5HMF)*(K_L_5HMF)*C)/(1+(K_L_5HMF)*C) ;
dqdt = kf*(qstar-q) ;

c=[1;1] ;
f = [De*dCdx;0] ;
s = [-u_s*dCdx-(1-epsilon)/epsilon*rho_liq*dqdt;dqdt] ;
end

%% initial condition
function u0 = icfun1(x)
u0 = [0;0] ;
end 

%% boundary condition
function [p1,q1,pr,qr] = bcfun1(x1,u1,xr,ur,t)

global C_5HMF

C0 = C_5HMF ;
p1 = [u1(1)-C0;0] ;
q1 = [0;1] ;
pr = [0;0] ;
qr = [1;1] ;
end

function[c,f,s] = pdefun2(x,t,u,dudx)

global qmax_furfural K_L_furfural kf De u_s Rep Dp epsilon rho_liq

C = u(1) ;
q = u(2) ;
dCdx = dudx(1) ;
dqdx = dudx(2) ;

qstar = ((qmax_furfural)*(K_L_furfural)*C)/(1+(K_L_furfural)*C) ;
dqdt = kf*(qstar-q) ;

c=[1;1] ;
f = [De*dCdx;0] ;
s = [-u_s*dCdx-(1-epsilon)/epsilon*rho_liq*dqdt;dqdt] ;
end

%% initial condition
function u0 = icfun2(x)
u0 = [0;0] ;
end 

%% boundary condition
function [p1,q1,pr,qr] = bcfun2(x1,u1,xr,ur,t)

global C_furfural

C0 = C_furfural ;
p1 = [u1(1)-C0;0] ;
q1 = [0;1] ;
pr = [0;0] ;
qr = [1;1] ;
end

function[c,f,s] = pdefun3(x,t,u,dudx)

global qmax_acetic K_L_acetic kf De u_s Rep Dp epsilon rho_liq

C = u(1) ;
q = u(2) ;
dCdx = dudx(1) ;
dqdx = dudx(2) ;

qstar = ((qmax_acetic)*(K_L_acetic)*C)/(1+(K_L_acetic)*C) ;
dqdt = kf*(qstar-q) ;

c=[1;1] ;
f = [De*dCdx;0] ;
s = [-u_s*dCdx-(1-epsilon)/epsilon*rho_liq*dqdt;dqdt] ;
end

%% initial condition
function u0 = icfun3(x)
u0 = [0;0] ;
end 

%% boundary condition
function [p1,q1,pr,qr] = bcfun3(x1,u1,xr,ur,t)

global C_acetic

C0 = C_acetic ;
p1 = [u1(1)-C0;0] ;
q1 = [0;1] ;
pr = [0;0] ;
qr = [1;1] ;
end

function[c,f,s] = pdefun4(x,t,u,dudx)

global qmax_phenol K_L_phenol kf De u_s Rep Dp epsilon rho_liq

C = u(1) ;
q = u(2) ;
dCdx = dudx(1) ;
dqdx = dudx(2) ;

qstar = ((qmax_phenol)*(K_L_phenol)*C)/(1+(K_L_phenol)*C) ;
dqdt = kf*(qstar-q) ;

c=[1;1] ;
f = [De*dCdx;0] ;
s = [-u_s*dCdx-(1-epsilon)/epsilon*rho_liq*dqdt;dqdt] ;
end

%% initial condition
function u0 = icfun4(x)
u0 = [0;0] ;
end 

%% boundary condition
function [p1,q1,pr,qr] = bcfun4(x1,u1,xr,ur,t)

global C_phenol

C0 = C_phenol ;
p1 = [u1(1)-C0;0] ;
q1 = [0;1] ;
pr = [0;0] ;
qr = [1;1] ;
end

function[c,f,s] = pdefunEB(x,t,u,dudx)

global qmax_phenol K_L_phenol kf De u_s Rep Dp epsilon rho_liq cp_AC cp_feed H_ads h alpha

C = u(1) ;
q = u(2) ;
Tf = u(3) ;
dCdx = dudx(1) ;
dqdx = dudx(2) ;
dTdx = dudx(3) ;

Ts = 298.15 ; % adsorbent temp Kelvin

qstar = ((qmax_phenol)*(K_L_phenol)*C)/(1+(K_L_phenol)*C) ;
dqdt = kf*(qstar-q) ;

c = ones(3,1) ;
f = [De*dCdx;
        0;
        0];
s = [-u_s*dCdx-(1-epsilon)/epsilon*rho_liq*dqdt;
        dqdt;
        -u_s*dTdx+(((h*alpha*(Tf-Ts)-(H_ads/cp_AC))+(H_ads/cp_feed))*(dqdt*((1-epsilon)/epsilon)))] ;
end

%% initial condition
function u0 = icfunEB(x)
u0 = [0;
        0;
        298.15] ;
end 

%% boundary condition
function [p1,q1,pr,qr] = bcfunEB(x1,u1,xr,ur,t)

global C_phenol

C0 = C_phenol ;
p1 = [u1(1)-C0;0;u1(3)-298.15] ;
q1 = [0;1;0] ;
pr = [0;0;0] ;
qr = [1;1;1] ;
end

[Torsten]

The error message changed - now ode15s seems to have problems integrating your system of ordinary differential equations (see above).

[Josh]

Yes, so the error message has changed since i implemented your boundary condition change... How would you go about fixing it ?

[Torsten]

As you can see above, the temperature explodes at the right endpoint (see above). Check the possible reasons (especially your sourceterm -u_s*dTdx+(((h*alpha*(Tf-Ts)-(H_ads/cp_AC))+(H_ads/cp_feed))*(dqdt*((1-epsilon)/epsilon))) )

h*alpha = 1.5e9 !!!