Having trouble with lsqcurvefit for a heat transfer experiment.

Hello all,

I am attempting to use lsqcurvefit in order to find values of h for a temperature distribution of a fin, I have 5 experimental temperature datapoints for each different part of the experiment. However, my values for h are nowhere near where they should be. I have been stuck on this for a few days and would appreciate any tips as I am unfamiliar with the lsqcurve fit function. Here is A2_TF.TXT

filename = 'A2_TF.TXT';
data = readmatrix(filename); % read the .csv file's data

%separating the main matrix into sets of data for each experiment
free_brass = data(152:838,1:5); 
forced_brass = data(1196:1882,1:5);
free_copper = data(152:838,6:10);
forced_copper = data(1196:1882,6:10);
free_steel = data(1196:1882,11:15);
forced_steel = data(152:838,11:15);
forced_alum = data(152:838,16:20);
free_alum = data(1196:1882,16:20);

%Only using the last row of data once the system has reach steady state, then converting to kelvin 
steadyfree_brass = 273 + free_brass(687,1:5);
steadyforced_brass = 273 + forced_brass(687,1:5);
steadyfree_copper = 273+ free_copper(687,1:5);
steadyforced_copper = 273+ forced_copper(687,1:5);
steadyfree_steel = 273+ free_steel(687,1:5);
steadyforced_steel = 273 + forced_steel(687,1:5);
steadyforced_alum = 273 + forced_alum(687,1:5);
steadyfree_alum = 273 + free_alum(687,1:5);

%Physical properties of the fin, here is where I begin to be unsure with
%difining h
h = 150;
D = 0.0127;
L = 0.305;
A = (pi.*D^2)/4;
P = pi.*D;
x = linspace(0,0.3048,5); %Length of the fin divided into 5 points, representing thermocouples 
Tinf = 295.2;

k_brass = 116;
k_copper = 385;
k_steel = 14.9;
k_aluminum = 167;

%Expression for m as a part of the temperature distribution of a fin.
m_brass = sqrt((h*P)/(k_brass*A));
m_copper = sqrt((h*P)/(k_copper*A));
m_steel = sqrt((h*P)/(k_steel*A));
m_aluminum = sqrt((h*P)/(k_aluminum*A));

%Theoretical values for temperature based on the temperature distribution
%formula, each fin is exposed to free and forced convection, that is the 2
%experiments.
T_free_brass = (steadyfree_brass(5)-Tinf)*(cosh(m_brass*(L-x))+(h/(m_brass*k_brass))*sinh(m_brass*(L-x)))/(cosh(m_brass*L)+(h/(m_brass*k_brass))*sinh(m_brass*L)) + Tinf

T_free_brass = 1×5

346.7094 306.2679 297.5822 295.7314 295.4055

T_free_copper = (steadyfree_copper(5)-Tinf)*(cosh(m_copper*(L-x))+(h/(m_copper*k_copper))*sinh(m_copper*(L-x)))/(cosh(m_copper*L)+(h/(m_copper*k_copper))*sinh(m_copper*L)) + Tinf

T_free_copper = 1×5

337.8449 313.6222 303.3240 299.1662 298.0062

T_free_steel = (steadyfree_steel(5)-Tinf)*(cosh(m_steel*(L-x))+(h/(m_steel*k_steel))*sinh(m_steel*(L-x)))/(cosh(m_steel*L)+(h/(m_steel*k_steel))*sinh(m_steel*L)) + Tinf

T_free_steel = 1×5

347.8527 295.9210 295.2099 295.2001 295.2000

T_free_aluminum = (steadyfree_alum(5)-Tinf)*(cosh(m_aluminum*(L-x))+(h/(m_aluminum*k_aluminum))*sinh(m_aluminum*(L-x)))/(cosh(m_aluminum*L)+(h/(m_aluminum*k_aluminum))*sinh(m_aluminum*L)) + Tinf

T_free_aluminum = 1×5

319.7510 302.0174 297.1016 295.7612 295.4758

T_forced_brass = (steadyforced_brass(5)-Tinf)*(cosh(m_brass*(L-x))+(h/(m_brass*k_brass))*sinh(m_brass*(L-x)))/(cosh(m_brass*L)+(h/(m_brass*k_brass))*sinh(m_brass*L)) + Tinf

T_forced_brass = 1×5

322.8562 301.1425 296.4790 295.4853 295.3104

T_forced_copper = (steadyforced_copper(5)-Tinf)*(cosh(m_copper*(L-x))+(h/(m_copper*k_copper))*sinh(m_copper*(L-x)))/(cosh(m_copper*L)+(h/(m_copper*k_copper))*sinh(m_copper*L)) + Tinf

T_forced_copper = 1×5

316.8976 304.5732 299.3335 297.2180 296.6278

T_forced_steel = (steadyforced_steel(5)-Tinf)*(cosh(m_steel*(L-x))+(h/(m_steel*k_steel))*sinh(m_steel*(L-x)))/(cosh(m_steel*L)+(h/(m_steel*k_steel))*sinh(m_steel*L)) + Tinf

T_forced_steel = 1×5

327.4872 295.6422 295.2061 295.2001 295.2000

T_forced_aluminum = (steadyforced_alum(5)-Tinf)*(cosh(m_aluminum*(L-x))+(h/(m_aluminum*k_aluminum))*sinh(m_aluminum*(L-x)))/(cosh(m_aluminum*L)+(h/(m_aluminum*k_aluminum))*sinh(m_aluminum*L)) + Tinf

T_forced_aluminum = 1×5

312.1585 299.9091 296.5135 295.5876 295.3905

%Here I try to use the lsqcurvefit function, but without any sucess because
%my values for h are around -2000 and -4000.
T_free_brass2 = @(h,x) (steadyfree_brass(5)-Tinf)*(cosh(m_brass*(L-x))+(h/(m_brass*k_brass))*sinh(m_brass*(L-x)))/(cosh(m_brass*L)+(h/(m_brass*k_brass))*sinh(m_brass*L)) + Tinf

T_free_brass2 = function_handle with value:
    @(h,x)(steadyfree_brass(5)-Tinf)*(cosh(m_brass*(L-x))+(h/(m_brass*k_brass))*sinh(m_brass*(L-x)))/(cosh(m_brass*L)+(h/(m_brass*k_brass))*sinh(m_brass*L))+Tinf

h_freebrass = lsqcurvefit(T_free_brass2,10,x,steadyfree_brass)

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

h_freebrass = -2.3318e+03

T_free_copper2 = @(h,x) (steadyfree_copper(5)-Tinf)*(cosh(m_copper*(L-x))+(h/(m_copper*k_copper))*sinh(m_copper*(L-x)))/(cosh(m_copper*L)+(h/(m_copper*k_copper))*sinh(m_copper*L)) + Tinf

T_free_copper2 = function_handle with value:
    @(h,x)(steadyfree_copper(5)-Tinf)*(cosh(m_copper*(L-x))+(h/(m_copper*k_copper))*sinh(m_copper*(L-x)))/(cosh(m_copper*L)+(h/(m_copper*k_copper))*sinh(m_copper*L))+Tinf

h_freecopper = lsqcurvefit(T_free_copper2,10,x,steadyfree_copper)

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

h_freecopper = -4.0230e+03

Answer 1

[Torsten]

The calls to lsqcurvefit must be

h_freebrass = lsqcurvefit(T_free_brass2,10,x,T_free_brass)
h_freecopper = lsqcurvefit(T_free_copper2,10,x,T_free_copper)

instead of

h_freebrass = lsqcurvefit(T_free_brass2,10,x,steadyfree_brass)
h_freecopper = lsqcurvefit(T_free_copper2,10,x,steadyfree_copper)

to reproduce the h value you are using.

Answer 2

[William Rose]

@Nate,

The plot below shows the measured data (red *), and the predicted temperature profile with h=150, 2000, and -4000. It seems that your data is not organized the way you expect, because the model, with "reasonable" h (i.e. h>0), does not look anything like the data. Perhaps the five data points are in reverse order. The plot also shows that there is no significant change in predicted temperature profile from h=150 to h=2000.

I like how you include a number of explanatory comments in your code. I recommend that yoou add mo4e comments. In particular, every time you specify a constant, say what the units are and what it represents. I have done that, with some guesses, below. What is h?

The code below shows free copper only, for clarity.

data = readmatrix('A2_TF.TXT'); % read the .csv file's data
%separating the main matrix into sets of data for each experiment
free_copper = data(152:838,6:10);
% Assume last row of data = steady state. Convert to kelvin. 
steadyfree_copper = 273+ free_copper(687,1:5);
TssCu5=steadyfree_copper(5);  % steady state copper temp at the tip (K)

% Fin properties. Here is where I begin to be unsure defining h
h = 150;            % what is h, and what units?
D = 0.0127;         % diameter (m?)
L = 0.305;          % length? (1 foot = 0.305 m?)
A = (pi.*D^2)/4;    % cross-section area (m^2?)
P = pi.*D;          % circumference (m?)
x = linspace(0,0.3048,5); %thermocouple locations (m) 
Tinf = 295.2;
k_copper = 385;
%Expression for m as a part of the temperature distribution of a fin.
m_copper = sqrt((h*P)/(k_copper*A));

%Theoretical values for temperature based on the temperature distribution
%formula, each fin is exposed to free and forced convection, that is the 2
%experiments.
T_free_copper = (TssCu5-Tinf)*(cosh(m_copper*(L-x))+(h/(m_copper*k_copper))*sinh(m_copper*(L-x)))/(cosh(m_copper*L)+(h/(m_copper*k_copper))*sinh(m_copper*L)) + Tinf;
%Here I try to use the lsqcurvefit function, but without any sucess because
%my values for h are around -2000 and -4000.
T_free_copper2 = @(h,x) (TssCu5-Tinf)*(cosh(m_copper*(L-x))+(h/(m_copper*k_copper))*sinh(m_copper*(L-x)))/(cosh(m_copper*L)+(h/(m_copper*k_copper))*sinh(m_copper*L)) + Tinf;
h_freecopper = lsqcurvefit(T_free_copper2,150,x,steadyfree_copper);

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

fprintf('Fitted h=%.2f.\n',h_freecopper)

Fitted h=-4022.96.

Plot x versus measured T; x versus Tpred (h=150); x vs Tpred (h=2000).

figure; plot(x,steadyfree_copper,'-r*',x,T_free_copper2(150,x),'-go',...
             x,T_free_copper2(2000,x),'-b+',x,T_free_copper2(-4000,x),'-mx'); 
grid on; xlabel('X (m)'); ylabel('Temp ({\circ}K)')
legend('Measured','Pred(h=150)','Pred(h=2000)','Pred(h=-4000)')

OK.