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

3 views (last 30 days)
Andrian Mirza
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
% 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
  2 Comments

Sign in to comment.

Answers (1)

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

Community Treasure Hunt

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

Start Hunting!