MATLAB Answers

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

조회 수: 1(최근 30일)
Andrian Mirza
Andrian Mirza 2021년 5월 11일
답변: 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
  댓글 수: 2
Andrian Mirza
Andrian Mirza 2021년 5월 12일
Please use n_tot instead of 1 in the loop!

댓글을 달려면 로그인하십시오.

답변(1개)

Chaitanya Mallela
Chaitanya Mallela 2021년 6월 23일
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!

Translated by