How to update a function with output from another function?

I am trying to solve for heat transfer coefficients so that I can use them to calculate gas and tank wall temperatures. Currently, I am using constant heat transfer coefficients to enable the calculation of the desired temperatures with ode45.

% Convective heat transfer coefficients
    h_in = 25; % Convective heat transfer coefficient inside the tank (W/(m^2·K))
    h_out = 6; % Convective heat transfer coefficient outside the tank (W/(m^2·K))

% Solve the coupled ODEs and PDE
    [t, y] = ode45(@(t, y) odes(t, y, R, C_p_gas, k_liner, k_CFRP, rho_liner, rho_CFRP, C_p_liner, C_p_CFRP, h_in, h_out, T_ambient, A_in, A_out, dx_liner, dx_CFRP, N_liner, N_CFRP, V, mdot_out), tspan, initial_conditions);
    
% Extract time histories for inside and outside wall temperatures
    T_wall_liner_interface = y(:,3); % Inside wall temperature (first point of the liner)
    T_wall_CFRP_interface = y(:,end); % Outside wall temperature (last point of the CFRP layer)

What I would like to do is to use the output temperatures from the ode45 to solve for h_in and h_out, which will in turn be used by ode45 to solve for new temperatures. The h_in and h_out calculations look like this:

    function [h_in, h_out]=heat_transfer_coefficients 
    Di = 0.230; % Internal diameter in meters
    Do = 0.279; % External diameter in meters

    T_gf = 266.17; %(K) - Type III hydrogen gas T at film temp
    rho_gf = 39.16; %(kg/m^3) - Type III hydrogen gas density @ T_gf

    Cp_gf = 0.01514*((T_gf/100)^4.789)+1002;
    mu_gf = (-4.127*10^-7)*((T_gf/100)^2)+((7.258*10^-6)*(T_gf/100))+(4.094*10^-7);
    lambda_gf = ((-3.77*10^-4)*(T_gf/100)^2)+((1.004*10^-2)*(T_gf/100))-(4.776*10^-4);
    beta_gf = 1/T_gf; %volumetric thermal expansion coefficient
    g = 9.81; %gravitational acceleration (m/s^2)

    Rai = (rho_gf^2*(Di^3)*beta_gf*g*(T_wall_liner_interface-T_gas)*Cp_gf)/(mu_gf*lambda_gf); %Rayleigh number at inner wall
    
    if Rai < 10^8
        Nu_i = 1.15*Rai.^0.22;
    else
        Nu_i = 0.14*Rai.^0.333;
    end

    ki = 168.9e3; %thermal conductivity inside tank at T_g, Type III (W/m*K)
    h_in = (Nu_i*ki)/Di;

    T_ambf = 279.43; %(K) - Type III
    rho_ambf = 1.263; %(kg/m^3) - density @ T_ambf
    c_wind = 0; %(m/s) - cross wind velocity

    Cp_ambf = 0.01514*((T_ambf/100)^4.789)+1002;
    mu_ambf = (-4.127*10^-7)*((T_ambf/100)^2)+((7.258*10^-6)*(T_ambf/100))+(4.094*10^-7);
    lambda_ambf = ((-3.77*10^-4)*(T_ambf/100)^2)+((1.004*10^-2)*(T_ambf/100))-(4.776*10^-4);
    beta_ambf = 1/T_ambf; %volumetric thermal expansion coefficient

    Re_Do = (rho_ambf*c_wind*Do)/(mu_ambf);
    Pr_ambf = (mu_ambf*Cp_ambf)/lambda_ambf;

    if Re_Do == 0
        Nu_cylforced = 0;
        Nu_sphforced = 0;
    else
        Nu_cylforced = 0.3+((0.62*(Re_Do^0.5)*(Pr_ambf^(1/3)))/(1+(0.4/((Pr_ambf))^(2/3)))^(1/4))...
        *(1+(Re_Do/282000)^5/8)^4/5;
        Nu_sphforced = 2 + ((Re_Do/4)+((3*10^-4)*Re_Do^1.6))^0.5;
    end

    Rao = (rho_ambf^2*(Do^3)*beta_ambf*g*(T_amb-T_wall_CFRP_interface)*Cp_ambf)/(mu_ambf*lambda_ambf); %Rayleigh number at inner wall
    
    Nu_oforced = ((Nu_cylforced*0.72)+(Nu_sphforced*0.24))/(0.72+0.24); 
    Nu_ofree = (0.6+((0.387*(Rao.^(1/6)))/(1+(0.559/Pr_ambf)^(9/16)).^(8/27))).^2; 
    Nu_o = ((Nu_ofree.^4)+(Nu_oforced^4)).^0.25;

    ko = 24.84e3; %thermal conductivity outside tank at T_wo, Type III (W/m*K)
    h_out = (Nu_o*ko)/Do;
    end

As can be seen in the above code, the h_in and h_out calculations require T_gas, T_wall_liner interface and T_wall_CFRP_interface, which currently are being calculated after h_in and h_out. Can somebody please show me a way I can incorporate the h_in and h_out calculations into the first section of code? For ease of explanation and solving, please feel free to use your own numbers for the constants in the ode45 solver line.

[Sam Chak]

But your function does not have inputs.

Note: Edited to avoid lengthy code.

[h_in, h_out] = heat_transfer_coefficients 

Unrecognized function or variable 'T_wall_liner_interface'.

Error in solution>heat_transfer_coefficients (line 16)
Rai = (rho_gf^2*(Di^3)*beta_gf*g*(T_wall_liner_interface-T_gas)*Cp_gf)/(mu_gf*lambda_gf); %Rayleigh number at inner wall

function [h_in, h_out]=heat_transfer_coefficients 
    ...
end

Answer 1

[Aquatris]

Here is one way iff the required variables are within the state vector 'y', then you can do:

%% change ode45 call
[t, y] = ode45(@(t, y) odes(t, y, R, C_p_gas, k_liner, k_CFRP,...
    rho_liner, rho_CFRP, C_p_liner, C_p_CFRP,...
    heat_transfer_coefficients(y),...% previously -> h_in, h_out,...
    T_ambient, A_in, A_out, dx_liner,...
    dx_CFRP, N_liner, N_CFRP, V, mdot_out), tspan, initial_conditions);

%% change odes function to allow input from heat_transfer_coefficient
function dy = odes(t, y, R, C_p_gas, k_liner, k_CFRP,...
    rho_liner, rho_CFRP, C_p_liner, C_p_CFRP,...
    h_in_h_out_vec,...% previously -> h_in, h_out,...
    T_ambient, A_in, A_out, dx_liner,...
    dx_CFRP, N_liner, N_CFRP, V, mdot_out)
    
    h_in = h_in_h_out_vec(1);
    h_out = h_in_h_out_vec(1);

    %% rest your code here %%
    ...
end

%% change heat_transfer_coefficients function to accept input arguments
function [h_in, h_out]=heat_transfer_coefficients(y)

    T_gas = y(1); % assuming it is the 1st element of your state vector, change accordingly
    T_wall_liner_interface = y(3);
    T_wall_CFRP_interface = y(end);
    %% rest of your code %%
    ...
end

[NewGuy]

Where is h_in_h_out_vec defined? Is this the output from the heat_transfer_coefficients function?

[NewGuy]

I've just tried running the code with your suggestions, but it's not working. The program is just running for a long time, without giving any results or error messages.

[Aquatris]

you need to provide your full code for us to diagnose. what is odes function? what are the values of the variables you use when you call odes within ode45?

My guess right now is it runs forever cause the ode is stiff or something that requires really small time steps to numerical integrate.

[NewGuy]

Please find attached my full code. I've made some changes to it since I last posted my question. Now, all the functions are in one script, and I've fixed some of the numerical errors I made.

Thank you for your help.

[Aquatris]

I was able to run your code and it gave the plots as well so I assume you dont have any issue anymore?

[Torsten]

@Thaw Nyi Nyi

Using the stiff integrator "ode15s" instead of "ode45" will be faster for the problem at hand.

[NewGuy]

@Aquatris I need to run the code for 3000 seconds but am unable to. In the code I uploaded, it was shorter than that since I was only testing it.

Thanks for your help and time. I really appreciate this since I’m still rather inexperienced with MATLAB.

[NewGuy]

@Torsten If I use ode15s, do I need to change anything in the code or can I directly replace ode45 with it, i.e. simply change ode45 wherever it appears to ode15s?

Thanks a lot for your help and patience. I’m still rather inexperienced with MATLAB so I really appreciate your advice.

[Torsten]

Yes, at the time where there is no hydrogen left in the tank, you'll get a problem with your equations. And this time is reached before 3000 s.

[NewGuy]

I've rechecked and determined I only need to run it for about 2750 s. Even then, the code just keeps running.

Can somebody please help me with this?

[Torsten]

Don't you read my responses ? Switch from ode45 to ode15s (no further changes in the call or elsewhere in the code are needed - just replace the names), and the solver will finish after about 2 seconds for tspan = [0, 2750].

[NewGuy]

Yes, I did read your comments. It was a misunderstanding on my part. I’m going to try out your suggestion.

Answer 2

[Torsten]

% Solve the coupled ODEs and PDE
[t, y] = ode45(@(t, y) odes(t, y, R, C_p_gas, k_liner, k_CFRP, rho_liner, rho_CFRP, C_p_liner, C_p_CFRP, T_ambient, A_in, A_out, dx_liner, dx_CFRP, N_liner, N_CFRP, V, mdot_out), tspan, initial_conditions);
 
function dy = odes(t, y, R, C_p_gas, k_liner, k_CFRP, rho_liner, rho_CFRP, C_p_liner, C_p_CFRP, T_ambient, A_in, A_out, dx_liner, dx_CFRP, N_liner, N_CFRP, V, mdot_out)
  T_gas = y(?)
  T_wall_liner_interface = y(3);
  T_wall_CFRP_interface = y(end);
  [h_in,h_out] = heat_transfer_coefficients(T_gas,T_wall_liner_interface,T_wall_CFRP_interface)
  ...
end

function [h_in, h_out]=heat_transfer_coefficients(T_gas,T_wall_liner_interface,T_wall_CFRP_interface) 
  ...
end

[Sam Chak]

Hi @Thaw Nyi Nyi

Let's test out @Torsten's concept of the solving the interdependent problem.

%% Main Script
tspan   = [0 30];   % duration of simulation
y0      = 1;        % initial value
[t, y]  = ode45(@(t, y) ode(t, y), tspan, y0);
plot(t, y), grid on
xlabel('t'), ylabel('y')

%% The hypothetical ODE (requires numerical solution of transcendental E)
function dy = ode(t, y)
    Esol= kepler(y);        % numerical solution of E (by calling kepler() function)
    dy  = - (y - Esol);     % ODE depends on the numerical solution of E
end

%% Kepler's equation (requires ODE solution y)
function Esol = kepler(M)
    e   = 0.5;                      % eccentricity of the orbit
    fun = @(E) M - (E - e*sin(E));  % cannot be solved for E algebraically. 
    E0  = 1;                        % initial guess (fixed, but can be made varying)
    opt = optimoptions('fsolve', 'Display', 'none');
    Esol= fsolve(fun, E0, opt);
end

[Sam Chak]

Hi @Thaw Nyi Nyi

The computation in the heat_transfer_coefficients() function can return values for h_in and h_out. However, you should be cautious, as certain input values can produce complex-valued solutions for h_in and h_out. Please ensure that the input values remain within the intended operating range.

%% Test 1
y             = [1, 1, 1, 1];
[h_in, h_out] = heat_transfer_coefficients(y)

h_in = 0

h_out = 3.2052e+04

%% Test 2
y             = [1, 2, 3, 4];
[h_in, h_out] = heat_transfer_coefficients(y)

h_in = 7.0213e+07 + 5.8086e+07i

h_out = 1.3315e+06 - 1.0367e+06i

function [h_in, h_out] = heat_transfer_coefficients(y)
    % extracting data from input
    T_wall_liner_interface  = y(1);
    T_gas                   = y(2);
    T_amb                   = y(3);
    T_wall_CFRP_interface   = y(4);

    % original code
    Di = 0.230; % Internal diameter in meters
    Do = 0.279; % External diameter in meters
    
    T_gf = 266.17; %(K) - Type III hydrogen gas T at film temp
    rho_gf = 39.16; %(kg/m^3) - Type III hydrogen gas density @ T_gf
    
    Cp_gf = 0.01514*((T_gf/100)^4.789)+1002;
    mu_gf = (-4.127*10^-7)*((T_gf/100)^2)+((7.258*10^-6)*(T_gf/100))+(4.094*10^-7);
    lambda_gf = ((-3.77*10^-4)*(T_gf/100)^2)+((1.004*10^-2)*(T_gf/100))-(4.776*10^-4);
    beta_gf = 1/T_gf; %volumetric thermal expansion coefficient
    g = 9.81; %gravitational acceleration (m/s^2)
    
    Rai = (rho_gf^2*(Di^3)*beta_gf*g*(T_wall_liner_interface-T_gas)*Cp_gf)/(mu_gf*lambda_gf); %Rayleigh number at inner wall
    
    if Rai < 10^8
        Nu_i = 1.15*Rai.^0.22;
    else
        Nu_i = 0.14*Rai.^0.333;
    end
    
    ki = 168.9e3; %thermal conductivity inside tank at T_g, Type III (W/m*K)
    h_in = (Nu_i*ki)/Di;
    
    T_ambf = 279.43; %(K) - Type III
    rho_ambf = 1.263; %(kg/m^3) - density @ T_ambf
    c_wind = 0; %(m/s) - cross wind velocity
    
    Cp_ambf = 0.01514*((T_ambf/100)^4.789)+1002;
    mu_ambf = (-4.127*10^-7)*((T_ambf/100)^2)+((7.258*10^-6)*(T_ambf/100))+(4.094*10^-7);
    lambda_ambf = ((-3.77*10^-4)*(T_ambf/100)^2)+((1.004*10^-2)*(T_ambf/100))-(4.776*10^-4);
    beta_ambf = 1/T_ambf; %volumetric thermal expansion coefficient
    
    Re_Do = (rho_ambf*c_wind*Do)/(mu_ambf);
    Pr_ambf = (mu_ambf*Cp_ambf)/lambda_ambf;
    
    if Re_Do == 0
        Nu_cylforced = 0;
        Nu_sphforced = 0;
    else
        Nu_cylforced = 0.3+((0.62*(Re_Do^0.5)*(Pr_ambf^(1/3)))/(1+(0.4/((Pr_ambf))^(2/3)))^(1/4))...
        *(1+(Re_Do/282000)^5/8)^4/5;
        Nu_sphforced = 2 + ((Re_Do/4)+((3*10^-4)*Re_Do^1.6))^0.5;
    end
    
    Rao = (rho_ambf^2*(Do^3)*beta_ambf*g*(T_amb-T_wall_CFRP_interface)*Cp_ambf)/(mu_ambf*lambda_ambf); %Rayleigh number at inner wall
    
    Nu_oforced = ((Nu_cylforced*0.72)+(Nu_sphforced*0.24))/(0.72+0.24); 
    Nu_ofree = (0.6+((0.387*(Rao.^(1/6)))/(1+(0.559/Pr_ambf)^(9/16)).^(8/27))).^2; 
    Nu_o = ((Nu_ofree.^4)+(Nu_oforced^4)).^0.25;
    
    ko = 24.84e3; %thermal conductivity outside tank at T_wo, Type III (W/m*K)
    h_out = (Nu_o*ko)/Do;
end