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

2 views (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
VBBV on 12 May 2021
How is it your program runs while loop? With given while loop condition
Andrian Mirza on 12 May 2021

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

### Categories

Find more on Ordinary Differential Equations in Help Center and File Exchange

### Community Treasure Hunt

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

Start Hunting!