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

Question

Andrian Mirza 2021-5-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/827815-dsolve-gives-an-error-when-trying-to-use-it-in-a-loop

回答： Chaitanya Mallela 2021-6-23

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