Errors using function and fsolve

Below is the code I have written and then I am getting the error messages. I am assuming my problem lies within the function or the fsolve. I have tried and tried and tried and cannot figure out what is wrong. Thank You for help!

I have tried so hard and would really be grateful for help! Thank You so much!

% *Set universal constants:*
sigma = 5.67e-8;    %  Stefan-Boltzmann constant (in W/m^2/K^4)
%
% *Set surrounding properties*
surr.T_flow = 295;     % surrounding flow temperature (in K)
surr.T_wall = 295;     % surrounding wall temperature (in K)
%
% *Set propeties for the computer chip*
chip.dy = 0.005;            % chip vertical thickness (in m)
chip.A_surf = 0.002;    % chip surface area in contact with heat sink (in m^2)
chip.k = 20;                  % chip thermal conductivity (in W/m/K)
%
% *Set propeties for the heat sink*
sink.eps = 0.7;         %  heat sink
sink.A_surf = 0.02; %  heat sink surface area  (in m^2) 
%%
% *Make vectors of convection coefficients and chip power*
 surr.h =[6,11,16,21,26]      %  heat transfer coefficient (in W/m^2*K)
n_h = length(surr.h);
 P_elec = [0:1:30]    %  chip power dissipation  (in W)
n_P = length(P_elec);
sink.T = zeros(n_P,n_h);  % allocate memory for sink temperatures (in K) 
for i_h=1:n_h
    for i_P=1:n_P
%%
% *Call the residual function to solve for sink temperature using f_zero*
   T_guess=300
      sink.T(i_h) = fsolve(@(T) Resfun(T,surr.h(i_h),sink,surr,sigma),T_guess);
    sink.T(i_P) = fsolve(@(T) Resfun(T,P_elec(i_P),sink,surr,sigma),T_guess);
      end
end
% *Residual function for finding sink temperature*
function res = Resfun(T,P_elec,sink,surr,sigma)%inputs  
 %  set up energy balance residual
  res=  P_elec - surr.h.*sink.A_surf.*(T-surr.T_flow)-sink.A_surf.*sigma.*sink.eps.*(T^4-surr.T_wall^4); % energy balance eq
end

[Brady Wayne Robinson]

There it is! Thank You

Answer 1

[Walter Roberson]

The function you are trying to find the zero of, must return a scalar. Your formula involves surr.h which is a vector, so your function is returning a vector.

Note that fzero which you used has this restriction, but fsolve does not have this restriction.

[Brady Wayne Robinson]

I tried using Fsolve and still got a boatload of errors

[Walter Roberson]

If you have the symbolic toolbox then

syms T
F = Resfun(T,P_elec(i_P),sink,surr,sigma)
S = solve(F(1),T);
vpa(subs(F(2:end), T, S))

You will see a matrix of non-zero values, too large in value to be just round off error. This tells you that the solution to the first of your equations cannot satisfy the other four equations. Therefore fsolve() is not going to be able to find a T value that solves all 5 equations simultaneously.

Your equations are hopeless in the form written.

[Walter Roberson]

Might I suggest to you that you do not want to try to find a T that satisfies all

 surr.h =[6,11,16,21,26]      %  heat transfer coefficient (in W/m^2*K)

simultaneously, and that instead you want to loop through solving your system for each one of those individually.

You seem to be giving that a try with

n_h = length(surr.h);
for i_h=1:n_h
...
        sink.T(i_h) = fsolve(@(T) Resfun(T,surr.h(i_h),sink,surr,sigma),T_guess);

so you are passing in the individual surr.h value as the second parameter. But inside Resfun you have

     res=  P_elec - surr.h.*sink.A_surf.*(T-surr.T_flow)-sink.A_surf.*sigma.*sink.eps.*(T^4-surr.T_wall^4); % energy balance eq

which uses all of surr.h

By the way,

sink.T(i_h) = fsolve(@(T) Resfun(T,surr.h(i_h),sink,surr,sigma),T_guess);

does not involve variables from the inner loop for i_P=1:n_P so there is no obvious point in doing inside the inner loop.

And

sink.T(i_P) = fsolve(@(T) Resfun(T,P_elec(i_P),sink,surr,sigma),T_guess);

is going to overwrite the same location at some point.

Perhaps you should be saving the results to a 2D array indexed at i_h and i_P . But you do need to resolve the question of why you are calling the function with two different inputs and writing to the same output variable.

Answer 2

[Ridwan Alam]

% *Set universal constants:*
global surr sink sigma

sigma = 5.67e-8;    %  Stefan-Boltzmann constant (in W/m^2/K^4)
%
% *Set surrounding properties*

surr.T_flow = 295;     % surrounding flow temperature (in K)
surr.T_wall = 295;     % surrounding wall temperature (in K)
%
% *Set propeties for the computer chip*
chip.dy = 0.005;            % chip vertical thickness (in m)
chip.A_surf = 0.002;    % chip surface area in contact with heat sink (in m^2)
chip.k = 20;                  % chip thermal conductivity (in W/m/K)
%
% *Set propeties for the heat sink*
sink.eps = 0.7;         %  heat sink
sink.A_surf = 0.02; %  heat sink surface area  (in m^2)

% *Make vectors of convection coefficients and chip power*
surr.h =[6,11,16,21,26];      %  heat transfer coefficient (in W/m^2*K)
n_h = length(surr.h);
P_elec = [0:1:30];    %  chip power dissipation  (in W)
n_P = length(P_elec);
sink.T = zeros(n_P,n_h);  % allocate memory for sink temperatures (in K)


for i_h=1:n_h
    for i_P=1:n_P

        % *Call the residual function to solve for sink temperature using f_zero*
        T_guess=300;
        temp_var = fsolve(@Resfun,[T_guess,surr.h(i_h)])
        sink.T(i_P,i_h) = temp_var(1);
        temp_var = fsolve(@Resfun,[T_guess,P_elec(i_P)])
        sink.T(i_P,i_h) = temp_var(1);
    end
end
% *Residual function for finding sink temperature*
function res = Resfun(alpha)%inputs  
    global surr sink sigma
    T = alpha(1);
    P_elec_i = alpha(2);
    %sink = alpha{3}; 
    %surr = alpha{4};
    %sigma = alpha{5};
    %  set up energy balance residual
    res=  P_elec_i - ...
        surr.h.*sink.A_surf.*(T-surr.T_flow) - ...
        sink.A_surf.*sigma.*sink.eps.*(T^4-surr.T_wall^4);
    % energy balance eq
end

Hope this helps.

[Walter Roberson]

global almost always causes more problems than it solves.

[Walter Roberson]

function res = Resfun(alpha)%inputs  
    global surr sink sigma
    T = alpha(1);
    P_elec_i = alpha(2);

Is the goal now to have Resfun vary T and P_elec_i to try to find a zero of the equation?

If so then calculate

           surr.h.*sink.A_surf.*(T-surr.T_flow) - ...
        sink.A_surf.*sigma.*sink.eps.*(T^4-surr.T_wall^4)

and set P_elec_i to that value, and you get a perfect 0 output.

Well, that is, you would, if you had fixed the problem that you are calculating for all surr.h values simultaneously, and I can guarantee you that you will not find a single [T, P_elect_i] pair that will solve all 5 equations simultaneously.