How to update a function with output from another function?

Question

NewGuy 2024-8-23，8:33

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2147439-how-to-update-a-function-with-output-from-another-function

编辑： NewGuy 2024-8-24，10:27

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Sam Chak 2024-8-23，10:22

编辑：Sam Chak 2024-8-23，10:23

在 MATLAB Online 中打开

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 2024-8-23，8:48

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2147439-how-to-update-a-function-with-output-from-another-function#answer_1504064

在 MATLAB Online 中打开

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

12 个评论
显示 10更早的评论隐藏 10更早的评论

Torsten 2024-8-23，21:33

编辑：Torsten 2024-8-23，22:06

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 2024-8-23，22:43

编辑：NewGuy 2024-8-24，10:27

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

请先登录，再进行评论。

Answer 2

Torsten 2024-8-23，10:32

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2147439-how-to-update-a-function-with-output-from-another-function#answer_1504199

编辑：Torsten 2024-8-23，10:33

在 MATLAB Online 中打开

% 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

2 个评论
显示无隐藏无

Sam Chak 2024-8-23，11:22

在 MATLAB Online 中打开

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 2024-8-23，11:37

在 MATLAB Online 中打开

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