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

Question

Josh 2024-4-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2109406-pdepe-solver-failure-spatial-discretization-partial-differential-derivative-temperature-modellin

编辑： Torsten 2024-4-20

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:

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Torsten 2024-4-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2109406-pdepe-solver-failure-spatial-discretization-partial-differential-derivative-temperature-modellin#answer_1444916

编辑：Torsten 2024-4-20

在 MATLAB Online 中打开

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.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Josh 2024-4-20

编辑：Torsten 2024-4-20

在 MATLAB Online 中打开

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

Josh 2024-4-20

编辑：Josh 2024-4-20

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

Torsten 2024-4-20

编辑：Torsten 2024-4-20

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

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

4 个评论
显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

4 个评论 显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论