Clear Filters
Clear Filters

How to update a function with output from another function?

3 views (last 30 days)
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 Comment
Sam Chak
Sam Chak on 23 Aug 2024
Edited: Sam Chak on 23 Aug 2024
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

Sign in to comment.

Answers (2)

Aquatris
Aquatris on 23 Aug 2024
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 Comments
Torsten
Torsten on 23 Aug 2024
Edited: Torsten on 23 Aug 2024
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
NewGuy on 23 Aug 2024
Edited: NewGuy on 24 Aug 2024
Yes, I did read your comments. It was a misunderstanding on my part. I’m going to try out your suggestion.

Sign in to comment.


Torsten
Torsten on 23 Aug 2024
Edited: Torsten on 23 Aug 2024
% 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 Comments
Sam Chak
Sam Chak on 23 Aug 2024
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
Sam Chak on 23 Aug 2024
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

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Tags

Products


Release

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!