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.

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:

%F Interp to stepsize 5  
x=0:5:l
F_fit=interp1(x_F,F,x)    

function ddx= odefun(~, x)  % x contains C 
% x(1:n)=C

%%% How to make C alonz z with ceratin value ?
% Interp F 
% C=rho./MW.*3.*(1-F_fit) 

    for j= 2:length(F_fit)   %since C= x(1:input.n)
    x(1:n,j)= rho./MW.*3.*(1-F_fit(j));   
   
    ddx(1)=MW.*(k0.*(rho).*exp(-Ea./(R.*x(n+1))).*x(1,j)...
        -A(1).*(rho).*exp(-Ea./(R.*x(n+1))).*(C_0-x(1,j)));
    
    
    ddx(2:n-1)=MW.*(k0.*(rho).*exp(-Ea./(R.*x(n+2:2*n-1))).*x(2:n-1,j)...
        -A(2:n-1).*(rho).*exp(-Ea./(R.*x(n+2:2*n-1))).*(C_0-x(2:n-1,j)));
    
    ddx(n)= MW.*(k0.*(rho).*exp(-Ea./(R.*x(2*n))).*x(n,j)...
        -A(n).*(rho).*exp(-Ea./(R.*x(2*n))).*(C_0-x(n,j)));
    
    end

ddx=ddx';
end

[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’.

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?

[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 made a vector x for the variables C and T, but it goes wrong....

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

Thanks in advance!!

[Torsten]

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

[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.

[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).

[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.