5.3 Nonlinear Curve Fitting

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+1  j, and h = hj+1  hj. 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/m2K. 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  eht) to the data.

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  eht and = teht

[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  eht)

{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.6710-8 W/m2K4

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)

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