Chapter 5
Curve Fitting
5.3 Nonlinear Curve Fitting
A sample consists of a layer of aluminum and a layer of a composite coating is tested in a vacuum chamber by measuring its temperature as a function of time. The behavior of the sample temperature has a nonlinear dependence on the emissivity of the sample and the heat transfer coefficient h between the coating and the vacuum air.
Figure 5.3-1. A sample enclosed within a testing chamber.
The unknown parameters h and may be obtained by fitting the model equation to experimental data as shown in Figure 5.3.2 where the curve represents the model equation and the circles represent the data.
Figure 5.3-2. Transient temperature of a typical sample.
To illustrate how this is done, first consider a portion of the graph in Figure 5.3-2 that is re-plotted in Figure 5.3-3. The relationship between the temperature Ti obtained from the model equation and the experimental value Ti,exp can be expressed generally as
Ti,exp = Ti (t; , h) + ei(5.3-2)
where ei is a random error that can be negative or positive. Ti is a function of the independent variable ti and the parameters h and . The random error is also called the residual, which is the difference between the calculated and measured values.
Figure 5.3-3. Relationship between the model equation and the data
Nonlinear regression is based on determining the values of the parameters that minimize the sum of the squares of the residuals called an objective function Fobj.
Fobj = = =(5.3-3)
Where N is the number of data points or measured temperatures in this case. The temperature from equation 5.3-1 can be expanded in a Taylor series around h and and curtailed after the first derivative.
Ti,j+1 = Ti,j + + h(5.3-4)
Where j is the guess and j+1 is the prediction, = j+1j, and h = hj+1hj. We have linearized the original model with respect to the parameters h and . Equation (5.3-4) can be substituted into Eq. (5.3-2) to yield
Ti,exp Ti,j = + h + ei(5.3-5a)
or in matrix form
{D} = [Zj]{A} + {E}(5.3-5b)
where [Zj] is the matrix of partial derivatives of the function (called the Jacobian matrix) evaluated at the guess j, the vector {D}contains the differences between the measure temperature and the calculated temperature at the guess j, the vector {A}contains the changes in the parameter values, and the vector {E} contains the residuals. It should be noted that as the final values of the parameters are obtained after the iterations vector {D} is the same as vector {E}.
[Zj] = , {D} = , {A} = , {E} =
We minimize the objective function
Fobj = = (5.3-3)
by taking its derivative with respect to each of the parameters and setting the resulting equation to zero.
= 2 = 0(5.3-6a)
= 2 = 0(5.3-6b)
This algorithm is Gauss-Newton method for minimizing the sum of the squares of the residuals between data and nonlinear functions. Equations (5.3-6a) and (5.3-6b) can be combined in a matrix form
[Zj]T{E} = 0(5.3-7)
where [Zj]T is the transpose of [Zj]. Let consider N = 3 so we can see the combination from (5.3-6a) and (5.3-6b) to (5.3-7).
= 0
Substitute {E} = {D} [Zj]{A} from Eq. (5.3-5b) into (5.3-7)
[Zj]T{{D} [Zj]{A}} = 0
or
[Zj]T[Zj]{A} = {[Zj]T{D}}(5.3-8)
The Jacobian matrix [Zj] may be evaluated numerically for the model equation (5.3-1).
(5.3-9a)
(5.3-9b)
Typically, can be chosen to be 0.01 and h can be chosen to be 0.01 W/m2K. Thus, the Gauss-Newton method consists of solving Eq. (5.3-8) for {A}, which can be employed to compute improved values for the parameters h and .
j+1 = j + (from {A})
hj+1 = hj + h (from {A})
This procedure is repeated until the solution converges that is until and h fall below an acceptable criterion. The Gauss-Newton method is a common algorithm that can be found in many numerical method texts. However, this description follows the notations and development by Chapra and Canale.
Example 5.3-1
Fit the function T(t; , h) = (1 eht) to the data.
t / 0.25 / 0.75 / 1.25 / 1.75 / 2.25T / 0.28 / 0.75 / 0.68 / 0.74 / 0.79
Use initial guesses of h = 1 and = 1 for the parameters.
Solution
The partial derivatives of the function with respect to the parameters h and are
= 1 eht and = teht
[Zj] = =
The matrix multiplied by its transpose results in
[Zj]T[Zj] =
[Zj]T[Zj] =
The vector {D} consists of the differences between the measurements T and the model predictions T(t; , h) = (1 eht)
{D} = =
The vector {D} is pre-multiplied by [Zj]T to give
[Zj]T{D} = =
The vector {A} can be calculated by using MATLAB statement dA=ZjTZj\ZjTD({A} = {[Zj]T[Zj] \[Zj]T{D}})
{A} =
The next guesses for the parameters h and are
h = 1 0.2715 = 0.7285
= 1 + 0.5019 = 1.5019
Table 5.3-1 lists the MATLAB program with the results of two iterations.
Table 5.3-1 ______
% Gauss-Newton method
%
t=[0.25 0.75 1.25 1.75 2.25]';
T=[0.28 0.57 0.68 0.74 0.79]';
e=1;h=1;
Tmodel='e*(1-exp(-h*t))';
dTde='1-exp(-h*t)';dTdh='e*t.*exp(-h*t)';
for i=1:2
Zj=[eval(dTde) eval(dTdh)];
ZjTZj=Zj'*Zj
D=T-eval(Tmodel)
ZjTD=Zj'*D
dA=ZjTZj\ZjTD
e=e+dA(1);
h=h+dA(2);
fprintf('Iteration #%g: e = %8.4f, h = %8.4f\n',i,e,h)
end
> Gauss
ZjTZj =
2.3194 0.9489
0.9489 0.4404
D =
0.0588
0.0424
-0.0335
-0.0862
-0.1046
ZjTD =
-0.1534
-0.0366
dA =
-0.2715
0.5019
Iteration #1: e = 0.7285, h = 1.5019
ZjTZj =
3.0660 0.4162
0.4162 0.0780
D =
0.0519
0.0777
0.0629
0.0641
0.0863
ZjTD =
0.2648
0.0397
dA =
0.0625
0.1758
Iteration #2: e = 0.7910, h = 1.6777
The parameters h and can also be calculated if the model equation is a differential equation. Consider the following nonlinear model that describes the temperature of the sample as a function of time.
(MaCp,a + McCp,c) = – Ah(T– Tair) – A( T4 – T4sur)(5.3-10)
where
Ma= mass of base aluminum
Mc= mass of composite coating
Cp,a= heat capacity of base aluminum
Cp,c= heat capacity of composite coating
A= surface area of coating
h= heat transfer coefficient between the air and the coating
= thermal emissivity of the coating
= Stefan-Boltzmann constant = 5.6710-8 W/m2K4
T= sample temperature assumed to be uniform over the aluminum and the coating
Tair= air temperature at the location below the sample
Tsur= surrounding temperature taken to be average of the bottom and side wall temperatures of the test chamber
The only unknown parameters in the differential equation (5.3-10) are h and . These parameters may be obtained by fitting the model equation to experimental data as shown in Figure 5.3-2 where the curve represents the model equation and the circles represent the data. It should be noted that once equation (5.3-10) is integrated, the calculated temperature Ti depends nonlinearly on its parameters h and , hence a nonlinear regression is required.
The iteration procedure is similar to the case when the model equation is given explicitly as a nonlinear equation. The Jacobian matrix however must be evaluated numerically using equations (5.3-9a, b)
(5.3-9a)
(5.3-9b)
The vector {D} consists of the differences between the measurements T and the model predictions T(t; , h) that must be integrated numerically. Table 5.3-2 lists a typical data for model equation (5.3-2) for which T = Tsample, Tsur = 0.5(Tbottom + Tside) where Tbottom is the temperature at the bottom surface and Tside is the temperature at the side surface of the test chamber. The temperatures are given at 1 minute interval.
Table 5.3-2 A typical data for model equation (5.3-10)
t(minute) / Tbottom / Tside / Tair / Tsample / t(minute) / Tbottom / Tside / Tair / Tsample1 / -80.9 / -74.8 / 50.6 / 65.5 / 31 / -100 / -95.5 / -67.6 / -51.4
2 / -80.9 / -67.4 / 52.7 / 62.5 / 32 / -99.3 / -94.9 / -68.8 / -53.4
3 / -81.3 / -71.4 / -31.9 / 55.6 / 33 / -102.1 / -97.3 / -70.3 / -55.2
4 / -80.7 / -72.4 / -47.6 / 46.4 / 34 / -103.2 / -98.4 / -71.5 / -57
5 / -79.8 / -71.9 / -28.7 / 38.3 / 35 / -103.2 / -98.4 / -72.4 / -58.6
6 / -78.9 / -71.4 / -22.8 / 32.4 / 36 / -102.9 / -97.6 / -73.3 / -60.4
7 / -77.7 / -70.6 / -19.7 / 26.7 / 37 / -101.8 / -96.9 / -73.9 / -62
8 / -76.4 / -70.5 / -19.7 / 21.3 / 38 / -100.5 / -96.4 / -74.6 / -63.4
9 / -75.3 / -70.3 / -21 / 16.5 / 39 / -99.4 / -95.7 / -75 / -64.9
10 / -74.8 / -70.1 / -22.4 / 11.4 / 40 / -98.5 / -94.9 / -75.3 / -66
11 / -75.1 / -70.3 / -24.2 / 7.3 / 41 / -97.6 / -94.4 / -75.7 / -67.6
12 / -73.9 / -70.1 / -26 / 3.2 / 42 / -96.7 / -93.9 / -75.9 / -68.5
13 / -73 / -69.7 / -27.8 / -0.4 / 43 / -96.9 / -94.4 / -76.6 / -69.9
14 / -72.1 / -69.4 / -29.4 / -3.9 / 44 / -98.2 / -97.3 / -77.5 / -70.8
15 / -76.4 / -73.7 / -31.9 / -7.1 / 45 / -99.3 / -99.1 / -78.9 / -71.7
16 / -83.6 / -77.1 / -34.5 / -10.3 / 46 / -100 / -100 / -79.8 / -72.6
17 / -87.7 / -83.4 / -37.3 / -13.8 / 47 / -100.3 / -100.3 / -80.4 / -73.7
18 / -89 / -87.7 / -42 / -16.8 / 48 / -100.3 / -100.2 / -81.1 / -74.8
19 / -89.9 / -89.7 / -44.7 / -19.9 / 49 / -100.2 / -99.1 / -81.4 / -75.5
20 / -93.3 / -92.2 / -46.7 / -22.9 / 50 / -99.6 / -97.8 / -81.6 / -76.4
21 / -99.1 / -94.4 / -50.3 / -25.8 / 51 / -99.3 / -97.3 / -81.8 / -77.3
22 / -101.4 / -94.8 / -53.4 / -28.5 / 52 / -98.5 / -96.6 / -82 / -78.2
23 / -107.4 / -99.4 / -55.5 / -31.6 / 53 / -97.8 / -96 / -82.3 / -78.9
24 / -107.5 / -100.9 / -58 / -34.3 / 54 / -97.5 / -96 / -82.3 / -79.5
25 / -107.2 / -101.1 / -60.2 / -37 / 55 / -97.3 / -96 / -82.7 / -80.4
26 / -106.3 / -100.2 / -62 / -39.7 / 56 / -96.7 / -95.5 / -82.7 / -80.9
27 / -105 / -99.1 / -63.8 / -42 / 57 / -96.6 / -95.3 / -82.9 / -81.4
28 / -103.8 / -97.6 / -65.1 / -44.7 / 58 / -96 / -95.1 / -82.9 / -81.7
29 / -102.3 / -96.7 / -66 / -47.1 / 59 / -95.7 / -94.9 / -82.9 / -82.3
30 / -101.2 / -95.8 / -67.2 / -49.4 / 60 / -95.7 / -94.8 / -83.2 / -82.9
61 / -97.3 / -95.8 / -83.8 / -83.4
Table 5.3-3 lists the MATLAB program to evaluate h and . The temperature data in Excel format can be read in to the program by the statement file= xlsread(fname). In this program
Jac = [dTdemis dTdhc] [Zj] =
Amat=Jac'*Jac [Zj]T[Zj]
TsamT {D} , h hc, emis
Bcol=Jac'* TsamT [Zj]T{D}
delA=Amat\Bcol {A} = {[Zj]T[Zj] \[Zj]T{D}}
The values obtained from the model equation are plotted in Figure 5.3-4 along with the experimental values. The model curve was obtained with the converged values of h and .
Table 5.3-3 ______
% Emissivity calculation
clear
% *********************** DATA BLOCK ************************
% Temperatures are in degree F, length in inch, density in kg/m3
% Cp in J/kg.K, hc in W/m2.K
%
sigma=5.67e-8;
fname='sample1';
file= xlsread(fname);
nb=12;
Tb=file(nb:61,1);
Ts=file(nb:61,2);
Tsur=(Tb+Ts)/2;
Tair=file(nb:61,3);
Tsam=file(nb:61,4);
%
den=1674.218;Cp=209.55; % Coating properties (adhesive and particles)
dena=2700;Cpa=900; %aluminum base properties
dt=60;
dt2=30;
% Th should be the initial temperature of the sample holder which is
% the same as the sample temperature
%
%
% *********************** END DATA BLOCK **********************
%
% Change temperature from F to K
Tsur=(Tsur+460)/1.8;Tair=(Tair+460)/1.8;Tsam=(Tsam+460)/1.8;
area=.001457;thick=0.000546;thicka=(1/16)*.0254;vol=area*thick;vola=area*thicka;
mCp=(vol*den*Cp)+(dena*vola*Cpa);con=sigma*area/mCp;conh=area/mCp;
%
Th=Tsam(1);
%
% Fourth Order Runge Kutta method
%
% Initial guess of h and emissivity
%
emis=.5;hc=2;
dhc=.01;demis=.01;
%
irow=length(Tair);tm=0:irow-1;
h=dt;hh=.5*h;h6=h/6;
Tcal=Tsam;
%
% Start curve fitting for hc and emissivity
%
for k=1:20
con1=emis*con;con2=hc*conh;T=Tsam(1);
for i=2:irow
Tsr=Tsur(i-1);
Ta=Tair(i-1);
k1=-con1*(T^4-Tsr^4)-con2*(T-Ta);
Tt=T+hh*k1;Ta=(Tair(i-1)+Tair(i))/2;
Tsr=(Tsur(i-1)+Tsur(i))/2;
k2=-con1*(Tt^4-Tsr^4)-con2*(Tt-Ta);
Tt=T+hh*k2;
k3=-con1*(Tt^4-Tsr^4)-con2*(Tt-Ta);
Tt=T+h*k3;Ta=Tair(i); Tsr=Tsur(i);
k4=-con1*(Tt^4-Tsr^4)-con2*(Tt-Ta);
T=T+h6*(k1+2*k2+2*k3+k4);
Tcal(i)=T;
end
TsamT=Tsam-Tcal;Tsav=Tcal;
% plot(tm,Tcal,tm,Tsam,'o')
% xlabel('Time(min.)');ylabel('Temperature(K)')
%
% Evaluate dT/dhc
%
hc=hc+dhc;
con2=hc*conh;T=Tsam(1);
for i=2:irow
Ta=Tair(i-1); Tsr=Tsur(i-1);
k1=-con1*(T^4-Tsr^4)-con2*(T-Ta);
Tt=T+hh*k1;Ta=(Tair(i-1)+Tair(i))/2;
Tsr=(Tsur(i-1)+Tsur(i))/2;
k2=-con1*(Tt^4-Tsr^4)-con2*(Tt-Ta);
Tt=T+hh*k2;
k3=-con1*(Tt^4-Tsr^4)-con2*(Tt-Ta);
Tt=T+h*k3;Ta=Tair(i);Tsr=Tsur(i);
k4=-con1*(Tt^4-Tsr^4)-con2*(Tt-Ta);
T=T+h6*(k1+2*k2+2*k3+k4);
Tcal(i)=T;
end
dTdhc=(Tcal-Tsav)/dhc;
%
% Evaluate dT/demis
%
hc=hc-dhc;emis=emis+demis;
con1=emis*con;con2=hc*conh;T=Tsam(1);
for i=2:irow
Ta=Tair(i-1);Tsr=Tsur(i-1);
k1=-con1*(T^4-Tsr^4)-con2*(T-Ta);
Tt=T+hh*k1;Ta=(Tair(i-1)+Tair(i))/2;
Tsr=(Tsur(i-1)+Tsur(i))/2;
k2=-con1*(Tt^4-Tsr^4)-con2*(Tt-Ta);
Tt=T+hh*k2;
k3=-con1*(Tt^4-Tsr^4)-con2*(Tt-Ta);
Tt=T+h*k3;Ta=Tair(i);Tsr=Tsur(i);
k4=-con1*(Tt^4-Tsr^4)-con2*(Tt-Ta);
T=T+h6*(k1+2*k2+2*k3+k4);
Tcal(i)=T;
end
dTdemis=(Tcal-Tsav)/demis;
Jac=[dTdemis dTdhc];
Amat=Jac'*Jac;Bcol=Jac'*TsamT;
delA=Amat\Bcol;
emis=emis-demis+delA(1);hc=hc+delA(2);
fprintf('e = %g, hc = %g\n',emis,hc)
if max(abs(delA))<.01, break, end
end
plot(tm,Tsav,tm,Tsam,'o')
xlabel('Time(min.)');ylabel('Temperature(K)')
e = 0.750943, hc = 2.54362
e = 0.295085, hc = 6.2165
e = 0.620875, hc = 4.09397
e = 0.620414, hc = 4.08901
Figure 5.3-4. Transient temperature of the sample given in Table 5.3-2.
1