Error using BVP4c (singular Jacobian encountered)

Question

Adarrsh Patil 2020-3-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/509684-error-using-bvp4c-singular-jacobian-encountered

回答： Devineni Aslesha 2020-4-23

Hi,

I am trying to get a concentration profile in my catalyst particle using the following code but I am getting an error of Singular Jacobian being encountered.

clear all

close all

clc

Pt = 85; %bar

P_CO = (1-0.05)/3*Pt; %bar

P_H2 = 2*(1-0.05)/3*Pt; %bar

P_H2O = 0*Pt;

P_CO2 = 0*Pt;

P_CH3OH = 0;

d_p = 3e-4; %diameter of particle in metres

xmesh = linspace(d_p/2,0,50);

solinit = bvpinit(xmesh, @guess);

sol = bvp4c(@bvpfcn, @bcfcn, solinit);

function dydx = bvpfcn(x,y)

P_CO = (1-0.05)/3*85; %bar

P_H2 = 2*(1-0.05)/3*85; %bar

P_CH3OH = 0;

T_rxn = 280; %degree Celsius

T = T_rxn + 273.15;

R = 8.314;

R1 = 8.314e-5; %bar.m^3.(mol.K)^-1

rho_cat = 1000; %kg/m3

D_CO = 1e-9; %m2/s

dydx = zeros(5,1);

dydx = [y(2)

(2.68e9*exp(-18400/T)*(K_calc(T)*y(1)*y(3)^2 - y(5)))*rho_cat*R1*T_rxn/(((1 + 0.069 * y(1) + 6.19e-8 * exp(6610/T) * y(3))^3)*D_CO*K_calc(T)) - 2*y(1)/x

y(3)

(2.68e9*exp(-18400/T)*(K_calc(T)*y(1)*y(3)^2 - y(5)))*rho_cat*R1*T_rxn*2/(((1 + 0.069 * y(1) + 6.19e-8 * exp(6610/T) * y(3))^3)*D_CO*K_calc(T)) - 2*y(3)/x

-y(1) - 2*y(3)];

function K_eq = K_calc(T)

A = 5139;

B = 12.621;

K_eq = 10 ^ (A/(T-B));

end

function res = bcfcn(ya,yb)

res = [ya(1) - (1-0.05)/3*85

yb(2)

ya(3) - (1-0.05)/3*85*2

yb(4)

ya(5)];