Grey box model of 1d heat diffusion in a rod (official mathworks example)

Hello,

I'm thinking about grey-box modelling of heat diffusion in a rod with the 1d heat equation.

There is a interesting example at mathworls side:

https://de.mathworks.com/help/ident/ug/estimate-continuous-time-grey-box-model-for-heat-diffusion.html?status=SUCCESS

I would adapt the example (image: original) and change it a little bit to include a a convection to air at the left end of the rod (image: adapted).

Here is the code for the original version with the insulation at one end:

function [A,B,C,D,K,x0] = heatd(kappa,htf,T,Ngrid,L,temp)
% ODE file parameterizing the heat diffusion model

% kappa (first parameter) - heat diffusion coefficient
% htf (second parameter) - heat transfer coefficient 
%                          at the far end of rod

% Auxiliary variables for computing state-space matrices:
% Ngrid: Number of points in the space-discretization
% L: Length of the rod
% temp: Initial room temperature (uniform)

% Compute space interval
deltaL = L/Ngrid;

% A matrix
A = zeros(Ngrid,Ngrid);
for kk = 2:Ngrid-1
  A(kk,kk-1) = 1;
  A(kk,kk) = -2;
  A(kk,kk+1) = 1;
end

% Boundary condition on insulated end
A(1,1) = -1; A(1,2) = 1;
A(Ngrid,Ngrid-1) = 1;
A(Ngrid,Ngrid) = -1;
A = A*kappa/deltaL/deltaL;

% B matrix
B = zeros(Ngrid,1);
B(Ngrid,1) = htf/deltaL;

% C matrix
C = zeros(1,Ngrid);
C(1,1) = 1;

% D matrix (fixed to zero)
D = 0;

% K matrix: fixed to zero
K = zeros(Ngrid,1);

% Initial states: fixed to room temperature
x0 = temp*ones(Ngrid,1);

This section should be altered, but I'm not sure how to do it to include a known htc for air at the right end.

% Boundary condition on insulated end
A(1,1) = -1; A(1,2) = 1;
A(Ngrid,Ngrid-1) = 1;
A(Ngrid,Ngrid) = -1;
A = A*kappa/deltaL/deltaL;

Then the idgrey object is constructed:

m = idgrey('heatd',{0.27 1},'c',{10,1,22});

Another thing I don't get is the structure of data. Unfortunately there is no explanation or example shown.

How can I implement my temperature data (for three points as two column table (t,T1) , (t,T2), (t,T3)) in the greyest?

me = greyest(data,m)

If I get it right, the parameters heat diffusion coefficient and heat transfer coefficient at the right end of rod are estimated, but I have to guess them as initial starting point? To do this the transient temperature data is required, unfortunately the example is incomplete.

Answer 1

[Torsten]

I think it will be simpler to couple "pdepe" for the solution of your PDE and "lsqnonlin" for optimizing your model parameters.

Maybe I'm wrong, but using "idgrey" for only 3 measurement points will force Ngrid = 3 which gives an insufficient resolution of the spatial temperature distribution.

Further, I find it hard to understand how to use "idgrey".

[Steffen B.]

Hello Thomas,

thanks for your advice. I guess you are right, in pdpe the issue is more tangible.

I changed the position of the points were measured data is available.

The liquid is now cooling the left side at x=0 m and the convction to air is at the right side located.

I used linear interpolation inside the bc function to implement a time dependent temperature values for the convection at the both ends of the rod (made of steel).

The code is looking like this:

clc,close all,clear all

L=1;
tsim=200;

nx=25;
dx=L/nx;
nt=200;
dt=tsim/nt;

x=linspace(0,L,nx);
t=linspace(0,tsim,nt);

Tdata=xlsread('TempData.xlsx');
save TempData.mat Tdata;

global time T_inf_fluid T_inf_air T1 T2 T3

time=Tdata(:,1);
T1=Tdata(:,3); % at x1=0; with fluid contact
T2=Tdata(:,5); % at x2=0.5; in the middle of the rod
T3=Tdata(:,7); % at x3=L; with air contact

T_inf_fluid=100;
T_inf_air=300;

m = 0;
sol = pdepe(m,@heatequation,@ic,@bc,x,t);

u = sol(:,:,:);

subplot(1,2,1)
surf(x,t,u,'EdgeColor','None')
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
view([150 25])

subplot(1,2,2)
plot(t,u(:,1))
hold on
plot(t,u(:,13))
hold on
plot(t,u(:,25))
hold on
xlabel('Time')
ylabel('Temperature u(0,t)')
title('Temperature change at xy')

function [c,f,s]=heatequation(x,t,u,dudx)
rho=7850;     % density
cp=420;       % heat capacity
lambda=45;    % conductivity
c=rho*cp;
f=lambda*dudx;
s=0;
end

function u0=ic(x)
Tini=300;      % initial temperature at t= 0 s
u0=Tini;
end

function [pl,ql,pr,qr]=bc(xl,ul,xr,ur,t)
global time T_inf_fluid T_inf_air T1 T2 T3
h_fluid=10000;
h_air=20;
pl=h_fluid*((interp1(time,T1,t,'linear','extrap'))-T_inf_fluid); % BC left, convection to a streaming fluid h_fluid=unknown
ql=1;                                                            %
pr=h_air*((interp1(time,T3,t,'linear','extrap'))-T_inf_air);     % BC right, convection to a streaming fluid h_air=known                                                           % 
qr=1;     
end

The result seems non-physical, a decereasing temperature field is expected due to the measured temperature values at the both ends of the rod. In the picture below is the calculated outcome (left side in the image) and the "measured temperature values at rods ends" as the input for interpolation. I also attached the corresponding xlsx file for the temperature interpolation.

The temperature is going up instead of decreasing, obviously I made a mistake in the bc function.

[Torsten]

I don't understand why you use your data within "pdepe".

From the optimizer, you will get passed values for lambda, h_fluid and h_air. With these values, the temperature profile over the rod is uniquely determined - thus also the temperatures in your three measurement points. By the difference between the simulated and measured temperatures, the optimizer will now try to adjust lambda, h_fluid and h_air as to make these differences as small as possible.

Thus your "pdepe" code would look like

L=1;
tsim=200;
nx=25;
dx=L/nx;
nt=200;
dt=tsim/nt;
x=linspace(0,L,nx);
t=linspace(0,tsim,nt);
Tini=300;     % initial temperature at t= 0 s
rho=7850;     % density
cp=420;       % heat capacity
lambda=45;    % conductivity
h_fluid=10000;
T_inf_fluid=100;
h_air=20;
T_inf_air=300;
m = 0;
sol = pdepe(m,@(x,t,u,dudx)heatequation(x,t,u,dudx,rho,cp,lambda),...
              @(x)ic(x,Tini),@(xl,ul,xr,ur,t)bc(xl,ul,xr,ur,t,h_fluid,T_inf_fluid,h_air,T_inf_air),x,t);
u = sol(:,:,:);
subplot(1,2,1)
surf(x,t,u,'EdgeColor','None')
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
view([150 25])
subplot(1,2,2)
plot(t,u(:,1))
hold on
plot(t,u(:,13))
hold on
plot(t,u(:,25))
hold on
xlabel('Time')
ylabel('Temperature u(0,t)')
title('Temperature change at xy')

●

function [c,f,s]=heatequation(x,t,u,dudx,rho,cp,lambda)
  c = rho*cp;
  f = lambda*dudx;
  s = 0;
end
function u0=ic(x,Tini)
  u0 = Tini;
end
function [pl,ql,pr,qr]=bc(xl,ul,xr,ur,t,h_fluid,T_inf_fluid,h_air,T_inf_air)
  pl = -h_fluid*(ul-T_inf_fluid); % BC left, convection to a streaming fluid h_fluid=unknown
  ql = 1;                                                            %
  pr = h_air*(ur-T_inf_air);     % BC right, convection to a streaming fluid h_air=known                                                           % 
  qr = 1;     
end

[Steffen B.]

Hello Thomas,

thank you very much for your helpful support.

In my model the temperature at the left and the right end of the rod is mandatory. Therefore I implemented the interpolation via interp1 (1D interpolation). With this, the temperature at the ends of the rod is corresponding to the data values. The convection coefficient h_air is known and therefore a constant. The material parameters are actual temperature dependent but I simplified it via known constants to keep it simple.

For now the only unknown value here is the h_fluid. The value should be optimized via lsqnonlin.

I'm struggling with the structure, it doesn't make sense to call lsqnonlin with the pedepe, pdepe has to calculate the temperature field for a guess h_fluid and then the field/ or specific points has to be compared with the corresponding experimental data Tdata respective T1, T2, T3.

clc,close all,clear all
L=1;
tsim=200;
nx=25;
dx=L/nx;
nt=200;
dt=tsim/nt;
x=linspace(0,L,nx);
t=linspace(0,tsim,nt);
Tdata=xlsread('TempData.xlsx');
save TempData.mat Tdata;
global time T_inf_fluid T_inf_air h_air rho cp lambda
p0=[1000];   % initial guess for h_fluid
lb=[10];     % lower value for h_fluid
ub=[100000]; % upper value for h_fluid
options = optimoptions('lsqnonlin','Algorithm','trust-region-reflective','Display','iter'); %  Guess of h_fluid with lsqnonlin fit
[p,resnorm,residual,exitflag,output]=lsqnonlin(@sol,p0,lb,ub,options) 
time=Tdata(:,1);
T1=Tdata(:,3); % at x1=x=0;
T2=Tdata(:,5); % at x2=x=0.5;
T3=Tdata(:,7); % at x3=x=L;
uexp=zeros(length(T1),3,1);
uexp=(T1,T2,T3);

Invalid expression. When calling a function or indexing a variable, use parentheses. Otherwise, check for mismatched delimiters.

TinterpT1=interp1(time,T1,t,'linear','extrap');
TinterpT3=interp1(time,T3,t,'linear','extrap');
T_inf_fluid=100;
T_inf_air=300;
h_air=20;
rho=7850;     % density
cp=420;       % heat capacity
lambda=45;    % conductivity
m= 0;
sol=pdepe(m,@heatequation,@ic,@bc(h_fluid),x,t);
u=sol(:,:,:);
ucomp1=u(:,1,1)-uexp(:,1); %  The structurte is still not right
ucomp2=u(:,13,1)-uexp(:,2);
ucomp3=u(:,25,1)-uexp(:,3);
f1=figure(1);
movegui(f1,[100 100]);
fontSize=20;
fontSize2=16;
f1.Position(3:4)=[1280 780];
subplot(1,2,1)
surf(x,t,u,'EdgeColor','None')
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
view([220 20])
LW=2;
subplot(1,2,2)
plot(t,u(:,1),"Color","blue","LineStyle","-","LineWidth",LW)
hold on
plot(t,u(:,13),"Color","black","LineStyle","-","LineWidth",LW)
hold on
plot(t,u(:,25),"Color","red","LineStyle","-","LineWidth",LW)
hold on
plot(time,T1,"Color","blue","LineStyle","--","LineWidth",LW)
hold on
plot(time,T2,"Color","black","LineStyle","--","LineWidth",LW)
hold on
plot(time,T3,"Color","red","LineStyle","--","LineWidth",LW)
hold on
xlabel('Time')
ylabel('Temperature u(0,t)')
title('Temperature change at xy')
function [c,f,s]=heatequation(x,t,u,dudx)
global rho cp lambda
c=rho*cp;
f=lambda*dudx;
s=0;
end
function u0=ic(x)
Tini=300;                               % initial temperature at t= 0 s
u0=Tini;
end
function [pl,ql,pr,qr]=bc(xl,ul,xr,ur,t,h_fluid)
global time T_inf_fluid T_inf_air h_air
pl=-h_fluid*(TinterpT1(t)-T_inf_fluid); % BC left, convection to a streaming fluid h_fluid=unknown
ql=1;                                                             %
pr=h_air*(TinterpT3(t)-T_inf_air);      % BC right, convection to a streaming fluid h_air=known                                                           % 
qr=1;                                                             %
end

The lsqnonlin should find the h_fluid with the least deviation concerning the three experimental points T1 (BC Left), T2 middle of the rod, and T3 (BC Right).

[Torsten]

You write the temperatures at both ends are mandatory, then you write you want to determine h_fluid with the least deviation concerning T1 and T3. Only one option can be correct.

If you want to prescribe the temperatures at both ends, you don't need h_fluid and h_air because you set temperature directly. In this case, no optimization is necessary and the way to validate your model is to compare T2 with the temperature obtained from "pdepe". The boundary condition function will set pl, ql, pr and qr as

pl = ul - TinterpT1(t)
ql = 0
pr = ur - TinterpT3(t)
qr = 0

If you want to get temperature profiles such that there might be a deviation from T1, T2 and T3, you need model constants that influence these temperatures and can be optimized. You suggest h_fluid, but I'm not sure if this single parameter will suffice. Setting

pl=-h_fluid*(TinterpT1(t)-T_inf_fluid); % BC left, convection to a streaming fluid h_fluid=unknown
ql=1;                                                             %
pr=h_air*(TinterpT3(t)-T_inf_air);      % BC right, convection to a streaming fluid h_air=known
qr=0;

makes no sense because for a mixed boundary condition, the temperature is a solution variable, namely ul and ur, and cannot be set to TinterpT1(t) and TinterpT2(t), respectively.

[Steffen B.]

Hello Thomas,

thanks for your helpfull support.

Yes you are right. Only the temperature or htc canbe defined at the boundary.

If I ran the code with the two temperatures at both ends i get different error messages here with the live script then at home.

clc,close all,clear all

L=1;
tsim=200;

nx=11;
dx=L/nx;
nt=200;
dt=tsim/nt;

x=linspace(0,L,nx);
t=linspace(0,tsim,nt);

Tdata=xlsread('TempData.xlsx');
save TempData.mat Tdata;

global TinterpT1 TinterpT3

time=Tdata(:,1);
T1=Tdata(:,3); % at x1=0;
T2=Tdata(:,5); % at x2=0.5;
T3=Tdata(:,7); % at x3=L;

TinterpT1=interp1(time,T1,t,'linear','extrap');
TinterpT3=interp1(time,T3,t,'linear','extrap');

m = 0;
sol=pdepe(m,@heatequation,@ic,@bc,x,t);

Array indices must be positive integers or logical values.

Error in solution>bc (line 72)
pl=ul-TinterpT1(t); % BC left, convection to a streaming fluid h_fluid=unknown

Error in pdepe (line 246)
[pL,qL,pR,qR] = feval(bc,xmesh(1),y0(:,1),xmesh(nx),y0(:,nx),t(1),varargin{:});

u = sol(:,:,:);

subplot(1,2,1)
surf(x,t,u,'EdgeColor','None')
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
view([130 20])

subplot(1,2,2)
plot(t,u(:,1),"Color","blue","LineStyle","-","LineWidth",1)
hold on
plot(t,u(:,5),"Color","black","LineStyle","-","LineWidth",1)
hold on
plot(t,u(:,11),"Color","red","LineStyle","-","LineWidth",1)
hold on
plot(time,T1,"Color","blue","LineStyle","--","LineWidth",1)
hold on
plot(time,T2,"Color","black","LineStyle","--","LineWidth",1)
hold on
plot(time,T3,"Color","red","LineStyle","--","LineWidth",1)
hold on
xlabel('Time')
ylabel('Temperature u(0,t)')
title('Temperature change at xy')

function [c,f,s]=heatequation(x,t,u,dudx)
rho=7850;     % density
cp=420;       % heat capacity
lambda=45;    % conductivity
c=rho*cp;
f=lambda*dudx;
s=0;
end

function u0=ic(x)
Tini=300;      % initial temperature at t= 0 s
u0=Tini;
end

function [pl,ql,pr,qr]=bc(xl,ul,xr,ur,t)
global TinterpT1 TinterpT3
pl=ul-TinterpT1(t);     % BC left, convection to a streaming fluid h_fluid=unknown
ql=0;                   %
pr=ur-TinterpT3(t);     % BC right, convection to a streaming fluid h_air=known                                                           % 
qr=0;     
end

At home I get:

Error in pdepe (line 250)

[pL,qL,pR,qR] = feval(bc,xmesh(1),y0(:,1),xmesh(nx),y0(:,nx),t(1),varargin{:});

Error in Test_only_Temperature (line 28)

sol=pdepe(m,@heatequation,@ic,@bc,x,t);

[Torsten]

Tdata=xlsread('TempData.xlsx');
time=Tdata(:,1);
T1=Tdata(:,3); % at x1=0;
T2=Tdata(:,5); % at x2=0.5;
T3=Tdata(:,7); % at x3=L;
TinterpT1 = @(t)interp1(time,T1,t);
TinterpT3 = @(t)interp1(time,T3,t);
L=1;
tsim=200;
nx=100;
dx=L/nx;
nt=200;
dt=tsim/nt;
x=linspace(0,L,nx);
t=linspace(1,tsim,nt);
Tini=300;     % initial temperature at t= 0 s
rho=7850;     % density
cp=420;       % heat capacity
lambda=45;    % conductivity
m = 0;
sol = pdepe(m,@(x,t,u,dudx)heatequation(x,t,u,dudx,rho,cp,lambda),...
              @(x)ic(x,Tini),@(xl,ul,xr,ur,t)bc(xl,ul,xr,ur,t,TinterpT1,TinterpT3),x,t);
u = sol(:,:,:);
plot(x,u(end,:))

●

function [c,f,s]=heatequation(x,t,u,dudx,rho,cp,lambda)
  c = rho*cp;
  f = lambda*dudx;
  s = 0;
end
function u0=ic(x,Tini)
  u0 = Tini;
end
function [pl,ql,pr,qr]=bc(xl,ul,xr,ur,t,TinterpT1,TinterpT3)
  pl = ul - TinterpT1(t); % BC left, convection to a streaming fluid h_fluid=unknown
  ql = 0;                                                            %
  pr = ur - TinterpT3(t); % BC right, convection to a streaming fluid h_air=known                                                           % 
  qr = 0;     
end

[Torsten]

As you can see from the results of the simulation with "pdepe", you will never get the T2 trajectory of your measurements when you use the partial differential equation as suggested by the model equation. The thermal mass of the rod is far too large to cool down within 200 seconds with the material data given. The plot above shows that after 200 seconds, the temperature in measurement point x2 still remains almost unchanged at Tini = 300 K.