dsolve gives an error when trying to use it in a loop

2 次查看(过去 30 天)
Can you run this code with L = 20000, n = 200?
It gives me an error. I am trying to discretize the function. So C_bulk0 should be substituded with C_bulk. It is something like this:
Can you suggest why this happens?
function [k_m, u, rho_G, D] = HgAds(L,n_tot)
%Constants
% Molar mass of Mercury
M_Hg = 200.59/1000; %kg/mol
% Avogadro's number
NA = 6.02*10^(23);
% Universal gas constant
R = 8.314; % J/(mol*K)
% Pipe length - L
% n - number of segments
% Temperature - T
% Pressure in - P_in Pascals
P_in = 143.2*10^(5); % Pa
% Pipe diameter
d_pipe = 0.7112; % meters
d_bulk = d_pipe - 4.52*10^(5);
% F_vol,in total
F_vol = 1.1586; % m^3/s
% Molar fraction of mercury - y0
% Methane - x0
% Molar Volume
%VG = R*T/P_in;
% Density of the gas
% M_CH4 = 16.042*10^(-3); % kg/mol
rho_G =0.68;
% Linear velocity
u = 2.2105;
% Viscosity muG -assuming only methane
%PcCH4=46.1e5; % Critical pressure
%TcCH4=190.6;
%N=0.0001778*(4.58*(T/TcCH4)-1.67)^0.625;
%mu_G =(4.6e-4*N*(M_CH4^(1/2)*PcCH4^(2/3))/(TcCH4^(1/6)))/1000;
mu_G = 2.15*10^(-5);
%Calculation of k_c, k_m
% Reynolds Number
Re = 1.98*10^(7);
% Schmidt Number
%D = (1.86*10^(-3)*T^1.5*(1/M_CH4 + 1/M_Hg)^0.5)/(P_in*3.364^2*1.578);
D = 9.17*10^(-8);
Sc = mu_G/rho_G/D;
Sh = 0.023*Re^0.8*Sc^0.33;
k_m = 0.0021;
SC = 0.95;
% DE1
syms C_bulk(t) C_stag(t) t
i = 1;
C_bulk0 = zeros(1,n_tot)+5000*10^(-12);
while i < 1
%Areas
A_stag = pi*d_bulk*L/n_tot;
A_pipe = pi*d_pipe*L/n_tot;
%Volumes
V_bulk = pi*d_bulk^2/4*L/n_tot;
V_pipe = pi*d_pipe^2/4*L/n_tot;
V_stag = V_pipe - V_bulk;
N_max = 1.2140*10^19;
Fug = 0.144992321;
s_0 = 1;
SSA = 1;
Z = 0.61;
v = 10^(17);
q_max = N_max*M_Hg/NA;
% ODES
ode1 = diff(C_bulk,t) == F_vol/V_bulk*(C_bulk0(i) - C_bulk) + k_m/V_bulk*A_stag*(C_stag - C_bulk);
%SC1 = (C_stag*Fug*Z*R*T1*NA/(M_Hg*sqrt(2*pi*M_Hg*R*T1))*s_0*(1-SC)/N_max - v*SC*exp(-1000*(151-28.82*SC)/(R*T1)));
ode2 = diff(C_stag,t) ==(k_m*A_stag*(C_bulk - C_stag) - SC*A_pipe*q_max*SSA)/V_stag;
odes = [ode1;ode2];
conds = [0,0];
[a,b] = vpa(dsolve(odes,conds));
i = i + 1;
C_bulk0(i) = a;
end
plot(a)
hold on
plot(b)
%odes = [ode1, ode2 ode3];
%[VF,Subs] = odeToVectorField(odes);
% odefcn = matlabFunction(VF, 'Vars',{t,Y});
%[t,y] = ode15s(odefcn,[0 1.167e+10],[0, 0, 0]);
%figure
%plot(t,y);
%grid;
%legend(string(Subs));
end

回答(1 个)

Chaitanya Mallela
Chaitanya Mallela 2021-6-23
Change the while loop condition in the code.

类别

Help CenterFile Exchange 中查找有关 Ordinary Differential Equations 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by