vector and for-loop issues in ode solver

Dear all, i am trying to solve a temperature-dependent kinetic expression in ode solver with concentration problem.

With presumption of C:

Q1: The code with the given temeparure works without putting in ode, but it doesn't work when it is for temperature simulation in ode solver.(Solved)

With the correct way to use C:

Q2: I've decided to solve mass and heat balance eqs. together. Since the concentration C of degas will be depentdent on F_degree along z direction, but ode solver will develop the z automatically --> I try to write a for loop for C with regarding F_degree in ode solver, but the code doesn't work at all....

Can anyone tell me where is the problem and how can i solve it? ? I'm not sure if making a loop for C in z direction(j) in ode is a correct way or not...

Any help or advice is greatly appreciated!!

Codes:

%Parameters
input.n=10;              % # in r 
input.Ea=63000;          % [J/mol]
input.R=8.314;           % [J/mol/K]
input.k0=505;            % [m^3/kg/s]
input.H=112500;          % [J/mol]
input.rho=360;           % [kg/m^3]
input.MW=0.3;            % [kg/mol]
input.T0=linspace(300,350,input.n)
l=0.8                    % [m] z len
p=2                      % [bar] 
input.r=linspace(0,0.2,input.n)  % radius=[m]

%Concentration
%%% Correct calculation %%%
%reaction:
% A-->A_X + 3*X

z_F=linspace(0,0.8,6)
F_degree=[0.04, 0.10, 0.25, 0.40, 0.65, 0.85]

C_0=input.rho./input.MW.*3.*(1-F_degree(1))  %[kg/m^3] % C_0: X gas still remained on substnace A at the start 
input.C_0=C_0
            
%%% plot 
for i=1:length(F_degree)
C(i)= input.rho./input.MW.*3.*(1-F_degree(i));         % C: X gas still remained on substance A   in the process
% remained 
end
C
plot(z_F,C,'r.',MarkerSize=10)
xlabel("z [m]"); ylabel("C [mol/m^3]")
title('C vs z')


%F_dgeree Interp to z stepsize 0.005  
z=0:0.005:l
F_degree_fit=interp1(z_F,F_degree,z)   % with real experimental data 
plot(z,F_degree_fit,'b*')
title("z vs F\_degree\_fit)")

% X gas equilibrium Rate     
% r = k0 *(rho)*exp(-Ea/(R*T))*C - k0/Keq *(rho)*exp(-Ea/(R*T))*(C_0-C)
% r = input.k0 *(input.rho)*exp(-input.Ea/(input.R*T))*input.C - input.k0./input.Keq *(input.rho)*exp(-input.Ea/(input.R*T))*(input.C_0-input.C)
% P.S. (C_0-C) means the released X gas 

x0= zeros(2.*input.n,1);
x0(1:input.n)=input.C_0               % C_0
x0(input.n+1:2*input.n)=input.T0      % Inlet temperature T0

[z, x] = ode45(@(z,x) odefun(z,x,input), 0:0.005:l, x0)   

[R,Z]=meshgrid(input.r,z)
surf(R,Z,x(:,input.n+1:2*input.n))

function ddz= odefun(~, x, input)  % x contains C and T 
% x(1:input.n)=C
% x(input.n+1:2*input.n)=T

p=2;                %[bar]
f_eq(1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(x(input.n+1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(x(input.n+1)-273.15)))));
f_eq(2:input.n-1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(x(input.n+2:2*input.n-1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(x(input.n+2:2*input.n-1)-273.15)))));
f_eq(input.n)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(x(2*input.n)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(x(2*input.n)-273.15)))));

Keq=f_eq./(1-f_eq) ;   % will be a vector of 1*input.n
Keq=Keq.';
%%% How to make C alonz z with ceratin value ????

% F_degree=[0.04, 0.10, 0.25, 0.40, 0.65, 0.85];
% Interp F_degree with stepsize: 0.005 --> F_degree_fit 
% C=input.rho./input.MW.*3.*(1-F_degree_fit) 

    for j= 2:length(F_degree_fit)   %since C= x(1:input.n)
    x(1:input.n,j)= input.rho./input.MW.*3.*(1-F_degree_fit(j));   
    
    
    
    dCdz(1)=input.MW.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+1))).*x(1,j)...
        -input.k0./Keq(1).*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+1))).*(input.C_0-x(1,j)));
    
    
    dCdz(2:input.n-1)=input.MW.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+2:2*input.n-1))).*x(2:input.n-1,j)...
        -input.k0./Keq(2:input.n-1).*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+2:2*input.n-1))).*(input.C_0-x(2:input.n-1,j)));
    
    dCdz(input.n)= input.MW.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(2*input.n))).*x(input.n,j)...
        -input.k0./Keq(input.n).*(input.rho).*exp(-input.Ea./(input.R.*x(2*input.n))).*(input.C_0-x(input.n,j)));
    
    
    dTdz(1) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+1))).*x(1,j)...
        -input.k0./Keq(1).*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+1))).*(input.C_0-x(1,j)));
    
    
    dTdz(2:input.n-1) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+2:2*input.n-1))).*x(2:input.n-1,j)...
        -input.k0./Keq(2:input.n-1).*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+2:2*input.n-1))).*(input.C_0-x(2:input.n-1,j)));
    
    
    dTdz(input.n) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(2*input.n))).*x(input.n,j)...
        -input.k0./Keq(input.n).*(input.rho).*exp(-input.Ea./(input.R.*x(2*input.n))).*(input.C_0-x(input.n,j)));
    
    end

ddz=[dCdz,dTdz]';
end

But i got this error:

Unrecognized function or variable 'F_degree_fit'.

Error in Question>odefun (line 64)
for j= 2:length(F_degree_fit)   %since C= x(1:input.n)

Error in Question (line 43)
[z, x] = ode45(@(z,x) odefun(z,x,input), 0:0.005:l, x0)   

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)
  odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

[Star Strider]

I do not see any assignment to ‘f_eq’ in the posted code.

It needs to exist in the calling script workspace before it is passed as an argument.

[uki71319]

Dear @Star Strider, for Q1, It needs to exist in the calling script workspace before it is passed as an argument.

Do you mean to put f_eq=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T-273.15))))) before the function or in odefun ??

But then f_eq has a variable T, which would be simulated in odefun, if i put it before odefun, then the error will say unrecozided T.

And even if I've put f_eq in odefun like this,

function dTdz= odefun(~,T,input,p,f_eq)

p=2
f_eq=1-(exp((0.08-7*10^-4*p).*(160+44.*log(p+3)-(T-273.15)))./(1+exp((0.08-7*10^-4.*p).*(160+44.*log(p+3)-(T-273.15)))))

f_eq(1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T(1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T(1)-273.15)))))
f_eq(2:input.n-1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T(2:input.n-1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T(2:input.n-1)-273.15)))))
f_eq(input.n)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T(input.n)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T(input.n)-273.15)))))

input.Keq=f_eq./(1-f_eq)    % will be a vector of 1*input.n

dTdz(1) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*T(1))).*input.C...
    -input.k0./input.Keq(1).*(input.rho).*exp(-input.Ea./(input.R.*T(1))).*(input.C_0-input.C));


dTdz(2:input.n-1) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*T(2:input.n-1))).*input.C...
    -input.k0./input.Keq(2:input.n-1).*(input.rho).*exp(-input.Ea./(input.R.*T(2:input.n-1))).*(input.C_0-input.C));


dTdz(input.n) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*T(input.n))).*input.C...
    -input.k0./input.Keq(input.n).*(input.rho).*exp(-input.Ea./(input.R.*T(input.n))).*(input.C_0-input.C));


dTdz=dTdz';
end

then i still get this error:

Unrecognized function or variable 'f_eq'.

Error in Question (line 32)
[z, T] = ode45(@(z,T) odefun(z,T,input,p,f_eq), 0:0.005:l, T0)     

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)
  odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

[Star Strider]

It appears that ‘f_eq’ needs to be an anonymous function.

It and all the arguments it needs will have to be passed to ‘odefun’.

[uki71319]

Dear @Star Strider, I then put the anonymous function for f_eq,

[z, T] = ode45(@(z,T) odefun(z,T,input), 0:0.005:l, T0)     

function dTdz= odefun(~,T,input)
p=2

f_eq=@(T) 1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T-273.15)))))

f_eq(1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T(1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T(1)-273.15)))))
f_eq(2:input.n-1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T(2:input.n-1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T(2:input.n-1)-273.15)))))
f_eq(input.n)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T(input.n)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T(input.n)-273.15)))))

Keq=f_eq./(1-f_eq) ;   % will be a vector of 1*input.n
Keq=Keq.'

dTdz(1) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*T(1))).*input.C...
    -input.k0./Keq(1).*(input.rho).*exp(-input.Ea./(input.R.*T(1))).*(input.C_0-input.C));


dTdz(2:input.n-1) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*T(2:input.n-1))).*input.C...
    -input.k0./Keq(2:input.n-1).*(input.rho).*exp(-input.Ea./(input.R.*T(2:input.n-1))).*(input.C_0-input.C));


dTdz(input.n) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*T(input.n))).*input.C...
    -input.k0./Keq(input.n).*(input.rho).*exp(-input.Ea./(input.R.*T(input.n))).*(input.C_0-input.C));


dTdz=dTdz';
end

but it turns out to be :

Conversion to function_handle from double is not possible.

Error in Question>odefun (line 40)
f_eq(1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T(1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T(1)-273.15)))))

Error in Question (line 34)
[z, T] = ode45(@(z,T) odefun(z,T,input), 0:0.005:l, T0)     

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)
  odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Answer 1

[Torsten]

%Parameters
input.n=10;              % # in r 
input.Ea=63000;          % [J/mol]
input.R=8.314;           % [J/mol/K]
input.k0=505;            % [m^3/kg/s]
input.H=112500;          % [J/mol]
input.rho=360;           % [kg/m^3]
input.MW=0.3;            % [kg/mol]
l=0.8                    % [m] z total length
p=2                      % [bar]

%Concentration
%%% Presumption_1st_try_in_ode, which is wrong 
F_degree=[0.04, 0.10, 0.25, 0.40, 0.65, 0.85]
% A--> B+3*gas 
input.C_0=input.rho./input.MW.*3.*(1-F_degree(1))    % [mol/m^3] the [C] of gas molecule still in A at the start
input.C=input.rho./input.MW.*3.*(1-F_degree(end))    % [mol/m^3] the [C] of gas molecule still in A

%%% Correct calculation for concentration %%%
z=linspace(0,0.8,6)
C_0=input.rho./input.MW.*3.*(1-F_degree(1))         % at the beginning 
C=input.rho./input.MW.*3.*(1-F_degree)
%%% plot 
for i=1:length(F_degree)
C(i)= input.rho./input.MW.*3.*(1-F_degree(i));   
end
C
plot(z,C,'r.',MarkerSize=10)
xlabel("z [m]"); ylabel("C [mol/m^3]")
title('C vs z')

%Rate
%r = k0 *(rho)*exp(-Ea/(R*T))*C - k0/Keq *(rho)*exp(-Ea/(R*T))*(C_0-C)
%r = input.k0 *(input.rho)*exp(-input.Ea/(input.R*T))*input.C - input.k0./input.Keq *(input.rho)*exp(-input.Ea/(input.R*T))*(input.C_0-input.C)

T0=linspace(300,350,input.n);
[z, T] = ode45(@(z,T) odefun(z,T,input), 0:0.005:l, T0)     

function dTdz= odefun(~,T,input)
p=2;
f_eq(1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T(1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T(1)-273.15)))));
f_eq(2:input.n-1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T(2:input.n-1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T(2:input.n-1)-273.15)))));
f_eq(input.n)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(T(input.n)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(T(input.n)-273.15)))));
Keq=f_eq./(1-f_eq);
Keq = Keq.';

dTdz(1) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*T(1))).*input.C...
    -input.k0./Keq(1).*(input.rho).*exp(-input.Ea./(input.R.*T(1))).*(input.C_0-input.C));
dTdz(2:input.n-1) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*T(2:input.n-1))).*input.C...
    -input.k0./Keq(2:input.n-1).*(input.rho).*exp(-input.Ea./(input.R.*T(2:input.n-1))).*(input.C_0-input.C));
dTdz(input.n) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*T(input.n))).*input.C...
    -input.k0./Keq(input.n).*(input.rho).*exp(-input.Ea./(input.R.*T(input.n))).*(input.C_0-input.C));
dTdz=dTdz';
end

[uki71319]

Dear @Torsten, sorry i just modify again. it works now!!

But what if in the future, i want to use f_eq as anoymous function in ode? How can i modify my code to work?

Also do you have any ideas for my Q2?

I was thinking that you told me to solve mass balance equation as well, but then the C concentration needs to solved or set at certain point on z direction..........

[Torsten]

But what if in the future, i want to use f_eq as anoymous function in ode? How can i modify my code to work?

Why ? I don't see any advantage to do so.

I was thinking that you told me to solve mass balance equation as well, but then the C concentration needs to solved or set at certain point on z direction

Yes, you get a set of discretized equations for C just like the ones you work with for T.

[uki71319]

Dear @Torsten,

Yes, you get a set of discretized equations for C just like the ones you work with for T.

I've edited my problem with the updated codes for mass and heat balance equatuons. I made a vector x for the variables C and T. As for how to demand C with certain F_degree on certain z point, i wrote a loop in ode, but it goes wrong....

Could you tell me how i can modify this issue...?

Thanks in advance!!

[uki71319]

Dear @Torsten, or should i write a nested function to define certain z point with certain C....??

I've interpolated the data of F_degree with the stepsize z for C, but the code doesn't work...

(P.S. C is the gas concentrain which is still remained in the A substance.

And thus C_0-C will be the concenterain of the released gas. )

function ddz= odefun(~, x, input)  % x contains C and T 
% x(1:input.n)=C
% x(input.n+1:2*input.n)=T

p=2;                %[bar]
f_eq(1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(x(input.n+1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(x(input.n+1)-273.15)))));
f_eq(2:input.n-1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(x(input.n+2:2*input.n-1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(x(input.n+2:2*input.n-1)-273.15)))));
f_eq(input.n)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(x(2*input.n)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(x(2*input.n)-273.15)))));

Keq=f_eq./(1-f_eq) ;   % will be a vector of 1*input.n
Keq=Keq.';
%%% How to make C alonz z with ceratin value ????
%%% or should i write a nested function????
z=0:0.005:l
F_degree_fit=interp1(z_F,F_degree,z)   % with real experimental data 

% F_degree=[0.04, 0.10, 0.25, 0.40, 0.65, 0.85];
% Interp F_degree with stepsize: 0.005 --> F_degree_fit 
% C=input.rho./input.MW.*3.*(1-F_degree_fit) 

    for j= 2:length(F_degree_fit)   %since C= x(1:input.n)
    x(1:input.n,j)= input.rho./input.MW.*3.*(1-F_degree_fit(j));   
    
    
    
    dCdz(1)=input.MW.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+1))).*x(1,j)...
        -input.k0./Keq(1).*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+1))).*(input.C_0-x(1,j)));
    
    
    dCdz(2:input.n-1)=input.MW.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+2:2*input.n-1))).*x(2:input.n-1,j)...
        -input.k0./Keq(2:input.n-1).*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+2:2*input.n-1))).*(input.C_0-x(2:input.n-1,j)));
    
    dCdz(input.n)= input.MW.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(2*input.n))).*x(input.n,j)...
        -input.k0./Keq(input.n).*(input.rho).*exp(-input.Ea./(input.R.*x(2*input.n))).*(input.C_0-x(input.n,j)));
    
    
    dTdz(1) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+1))).*x(1,j)...
        -input.k0./Keq(1).*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+1))).*(input.C_0-x(1,j)));
    
    
    dTdz(2:input.n-1) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+2:2*input.n-1))).*x(2:input.n-1,j)...
        -input.k0./Keq(2:input.n-1).*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+2:2*input.n-1))).*(input.C_0-x(2:input.n-1,j)));
    
    
    dTdz(input.n) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(2*input.n))).*x(input.n,j)...
        -input.k0./Keq(input.n).*(input.rho).*exp(-input.Ea./(input.R.*x(2*input.n))).*(input.C_0-x(input.n,j)));
    
    end

ddz=[dCdz,dTdz]';
end

[Torsten]

I cannot decipher which PDEs you are solving for T and C.

[uki71319]

Dear @Torsten,

Sorry for the cunfususion, they are besed on this equations:

But i've simplied the eqautions, thus i didn't put other parameters.

My question is since C is depedent on F_degree, which is the experimental data. I was thinking, instead wrting a for loop in odefun, is it possible to solve this problem by interpolating C from F_degree data based on z stepsize to get C_fit first, then wrting a nested function where x(C,T) to demand x(:,1)=C_fit?

Then T will also need to be solved according to the C_fit. But i got another issue again.

I try this way to code in this way:

%Parameters
input.n=10;              % # in r 
input.Ea=63000;          % [J/mol]
input.R=8.314;           % [J/mol/K]
input.k0=505;            % [m^3/kg/s]
input.H=112500;          % [J/mol]
input.rho=360;           % [kg/m^3]
input.MW=0.3;            % [kg/mol]
input.T0=linspace(300,350,input.n)
l=0.8                    % [m] z len
p=2                      % [bar] 
input.r=linspace(0,0.2,input.n)  % radius=[m]

z_C=linspace(0,0.8,6);
F_degree=[0.04, 0.10, 0.25, 0.40, 0.65, 0.85];
C_0=input.rho./input.MW.*3.*(1-F_degree(1))         % start point C_0 
input.C_0=C_0

z=0:0.005:l;
C=input.rho./input.MW.*3.*(1-F_degree); % with real experimental data 
C_fit=interp1(z_C,C,z);

x0= zeros(2.*input.n,1);
x0(1:input.n)=input.C_0           ;    % C_0
x0(input.n+1:2*input.n)=input.T0  ;    % Inlet temperature


[z, x] = ode45(@(z,x) odefun(z,x,input), 0:0.005:l, x0)  ; 
C_fit=x(:,1) ;
T=x(:,2)
 
%plot
[R,Z]=meshgrid(input.r,z)
surf(R,Z,T')    

    function ddz= odefun(~, x, input)
%         x(:,1)=C
%         x(:,2)=T
        p=2;                %[bar]
        f_eq(1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(x(input.n+1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(x(input.n+1)-273.15)))));
        f_eq(2:input.n-1)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(x(input.n+2:2*input.n-1)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(x(input.n+2:2*input.n-1)-273.15)))));
        f_eq(input.n)=1-(exp((0.08-6*10^-4*p).*(160+44.*log(p+3)-(x(2*input.n)-273.15)))./(1+exp((0.08-6*10^-4.*p).*(160+44.*log(p+3)-(x(2*input.n)-273.15)))));

        Keq=f_eq./(1-f_eq) ;   % will be a vector of 1*input.n
        Keq=Keq.';
        % rate
        %r = k0 *(rho)*exp(-Ea/(R*T))*C - k0/Keq *(rho)*exp(-Ea/(R*T))*(C_0-C) 
        % --> r = input.k0 *(input.rho)*exp(-input.Ea/(input.R*T))*input.C - input.k0./input.Keq *(input.rho)*exp(-input.Ea/(input.R*T))*(input.C_0-input.C)

        dTdz(1) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+1))).*x(1)...
            -input.k0./Keq(1).*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+1))).*(input.C_0-x(1)));


        dTdz(2:input.n-1) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+2:2*input.n-1))).*x(2:input.n-1)...
            -input.k0./Keq(2:input.n-1).*(input.rho).*exp(-input.Ea./(input.R.*x(input.n+2:2*input.n-1))).*(input.C_0-x(2:input.n-1)));


        dTdz(input.n) = -input.H.*(input.k0.*(input.rho).*exp(-input.Ea./(input.R.*x(2*input.n))).*x(input.n)...
            -input.k0./Keq(input.n).*(input.rho).*exp(-input.Ea./(input.R.*x(2*input.n))).*(input.C_0-x(input.n)));

        

        ddz=dTdz';     

    end

Then i got this error:

    Error using odearguments
@(Z,X)ODEFUN(Z,X,INPUT) returns a vector of length 10, but the length of initial conditions vector is 20. The vector returned by @(Z,X)ODEFUN(Z,X,INPUT) and the initial
conditions vector must have the same number of elements.

Error in ode45 (line 107)
  odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
  

Thus, i wonder if i should still stick to solve both T and C balance PDE toghether?

But then again, since i've already defined the C_fit value along z that i want, so the main issue is how to make T to be developed based on C_fit.....?

Thanks in advance!

[Torsten]

I described the way how I would determine the amount of developing gas in mol/(m^3*s) for each cell in z- and r- direction. This value can be used in the energy equation for each cell to compute the source term.

I cannot understand your C computations in the new code.

[uki71319]

Dear @Torsten,

I described the way how I would determine the amount of developing gas in mol/(m^3*s) for each cell in z- and r- direction. This value can be used in the energy equation for each cell to compute the source term.

doesn't that mean i need to give up using odesolver? Or is it still using the MOL with ode solver to solve it?

From what you've descibed last time:

You know the molar flux entering (ndot_in) and leaving (ndot_out) the reactor and the residence time T in the reactor.

Thus (ndot_in - ndot_out)*T [moles] of your substance degas.

You know the free volume V of the reactor.

Thus (ndot_in - ndot_out)*T/V [moles/m^3] degas in total.

Now using your F, you can distribute this total amount (ndot_in - ndot_out)*T/V to the different z-zones of the reactor.

For z-Zone i : (ndot_in - ndot_out)*T/V * deltaFi.

where I assume that the deltaFi sum to 1 (maybe in your case, you have to set ndot_out = 0).

Now you can distribute this total amount for z-zone i to the r-zones in z-Zone i:

For r-Zone j in z-Zone i : (ndot_in - ndot_out)*T/V * deltaFi * pi*(r_(j+1)^2-r_j^2) / (pi*R^2)

Now multiply this expression with -delta_r H and you'll arrive at an approximate value for the heat released in cell (i,j).

--> I thought this is the way for not using kinetic expression; With kinetic expression term, then i can easily use MOL to solve my equations. Am i wrong? sorry i get confused by your way...

Cuz if i use the for loop to directly solve PDEs for z zone first, then r zone --> Wouldn't it be difficult to apply for the radius boundary conditions.?....Thus, i wonder are there any possible ways to mantain my MOL codes and adjusting my T from C by just creating a nested function or a for loop for only for z?

My new idea is: since from my experimental data, C is already known for certain points along z based on F_degree, thus i just need to interpolate the C data into C_fit along z with z_stepsize. And thus i think i just need to deal with heat balance eq.

So i make the x variables contain C and T, thus with the solution in ODE: [z,x]

x(:,1:input.n)=C --> which i want to make it as C_fit.

x(:,input.n+1:2*input)=T

But then i am not sure how to modify the dTdz with C in reaction rate term, since:

%r = k0 *(rho)*exp(-Ea/(R*T))*C - k0/Keq *(rho)*exp(-Ea/(R*T))*(C_0-C) 

I am really stuck in this point....

[Torsten]

doesn't that mean i need to give up using odesolver? Or is it still using the MOL with ode solver to solve it?

It means that you use your measurement data to solve the energy equation using the MOL approach.

--> I thought this is the way for not using kinetic expression; With kinetic expression term, then i can easily use MOL to solve my equations. Am i wrong? sorry i get confused by your way...

So you know the kinetic parameters for the degassing reaction ? And why then do you need F_degree in your model ? If you have the kinetic parameters, you can calculate F_degree by solving simultaneously the equation for C and T.

My new idea is: since from my experimental data, C is already known for certain points along z based on F_degree, thus i just need to interpolate the C data into C_fit along z with z_stepsize. And thus i think i just need to deal with heat balance eq.

But this is exactly what I suggested: deduce the reaction term in mol/(m^3*s) from your measurement data prior to your simulation and only solve the energy equation with these values fixed. But why do you need K_eq / f_eq for this ?

As far as I see, you have two alternatives:

If you have the kinetic parameters for a general degassing kinetic, you don't need F_degree. You can solve C and T simultaneously and get F_degree from the simulation (assuming that the other model parameters are known and don't need to be determined).

If you don't have the kinetic parameters, then - prior to your simulation - you can estimate the molar fluxes for z (and assuming uniform distribution over r, also for r) from F_degree. Having these terms, you can include them in the energy equation and solve it alone (without using the equation for C).

And with respect to the mentionned aim of your simulation:

I find it quite strange that you want to fit the two parameters you mentionned (thermal conductivity and heat transfer coefficient) by some kind of simulation. There are well-suited experimental designs (especially to measure thermal conductivity).

[uki71319]

Dear @Torsten, thank you very much!!!!!!! I got your points now!

1. If you have the kinetic parameters for a general degassing kinetic, you don't need F_degree. You can solve C and T simultaneously and get F_degree from the simulation (assuming that the other model parameters are known and don't need to be determined).

Yes, then i computed this altrady. Then i just need to use parameter estimation (lsqcurve) to fit the experimental temperature and get the right thermal conductivity and wall heat coefficient.

2. If you don't have the kinetic parameters, then - prior to your simulation - you can estimate the molar fluxes for z (and assuming uniform distribution over r, also for r) from F_degree. Having these terms, you can include them in the energy equation and solve it alone (without using the equation for C).

This second way i need to think about how should i code it properly. So as you said, i will need to reduce my heat equation with F_degree to n mol/(m^3*s) (the molar fluxes for z) first. But prior to my simulation --> That means i need to write a for-loop? If i got time, i will think more about the right coding way to do it!

[Torsten]

But prior to my simulation --> That means i need to write a for-loop?

I think pencil and paper is more suited regarding the number of elements in F_degree.

[uki71319]

Dear @Torsten, thank you for your help.

Regarding to some questions:

But why do you need K_eq / f_eq for this?

Sorry i forgot to mention that this is a reversible dgassing reaction.....

But then I also realize i got a bug fot my rate ri on mass equation: for dCdz= (D/v)*(1/r)*d(r*dC/dr)dr + (1/v)*ri.

If it's a reversible equation, then C rate should be written as:

%r = -k0 *(rho)*exp(-Ea/(R*T))*C + k0/Keq *(rho)*exp(-Ea/(R*T))*(C_0-C) 

instead of

%r = k0 *(rho)*exp(-Ea/(R*T))*C - k0/Keq *(rho)*exp(-Ea/(R*T))*(C_0-C) 

Am i right? Also the Mw moleculor weight should be deleted.

[uki71319]

Dear @Torsten, so following your second suggestion,

You know the molar flux entering (ndot_in) and leaving (ndot_out) the reactor and the residence time T in the reactor.

Thus (ndot_in - ndot_out)*T [moles] of your substance degas.

You know the free volume V of the reactor.

Thus (ndot_in - ndot_out)*T/V [moles/m^3] degas in total

--> since i could directly use the experimental C [mol/m^3] from degree of F along z, so these sentesnces (ndot_in - ndot_out)*T/V [moles/m^3] turns out to be : CA_in - CA_out, [mol/m^3] right?

But then,

Now using your F, you can distribute this total amount (ndot_in - ndot_out)*T/V to the different z-zones of the reactor.

For z-Zone i : (ndot_in - ndot_out)*T/V * deltaFi.

where I assume that the deltaFi sum to 1 (maybe in your case, you have to set ndot_out = 0).

Now you can distribute this total amount for z-zone i to the r-zones in z-Zone i:

For r-Zone j in z-Zone i : (ndot_in - ndot_out)*T/V * deltaFi * pi*(r_(j+1)^2-r_j^2) / (pi*R^2)

-> How come these above make it to the final reaction unit: [mol/(m^3*s)] ? Since F_degree is no unit.

Now multiply this expression with -delta_r H and you'll arrive at an approximate value for the heat released in cell (i,j)

--> and later if I keep going to multiply (CA_in - CA_out) with (-delta_r H),

then the final unit become: [mol/(m^3)]*[J/mol] = [J/m^3] , instead of what we want: [J/(m^3*s)]...

Thank you for making these detail instructions before, but when i follow it throgh, i found it still confusing. Could you please tell me more about it?