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

1 view (last 30 days)
Andrian Mirza on 11 May 2021
Answered: Chaitanya Mallela on 23 Jun 2021
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
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
##### 2 CommentsShow 1 older commentHide 1 older comment
Andrian Mirza on 12 May 2021

Chaitanya Mallela on 23 Jun 2021
Change the while loop condition in the code.