%   Script KENNWRT.M
%
%  Kennwertermittlung
%  Beispiel: Waermetauscher
%  Eingangsgröße: T20 (Eintrittstemperatur)
%  Ausgangsgröße: T1 (Austrittstemperatur)
%
% J. Lunze
% 2.5.2003            30.4.2002
% für 2. Auflage: 23.6.2007
% für 3. Auflage: 19.1.2012
% für 4. Auflage: 18.12.2015
% für 5. Auflage: 12.11.2019
%
%
clear
close all
%                  Parameter
%
w1=2;               %  kg/min
w2=3;
V=2;                %  m^3
T10=85;         % Grad C
T20=20;         % Grad C
c1=1;               % cal/g/Grad
c2=c1;
m1=2;               % kg
m2=1;               % kg
A1=1*0.02*pi;   % m^2  Rohr von 1 m Länge, 2 cm Durchmesser
k=100;          % cal/min/m^2/Grad

n=2;
A=[-(k*A1+w1*c1)/m1/c1  k*A1/m1/c1;
   k*A1/m2/c2           -(k*A1+w2*c2)/m2/c2];
B=[0; w2/m2];
C=[1 0];
D=0;
Waermetauscher=ss(A, B, C, D);
%
%   Statisches Verhalten des Wärmetauschers
ks=dcgain(Waermetauscher)


%%   Übergangsfunktion des Waermetauschers
%
%
figure(1)
subplot(2,1,1)
[Y, Ty]=step(Waermetauscher);
plot(Ty, Y, 'b');
%axis([-0.005 0.035 -2 35]);
yticks([0 0.2 0.4 0.6]);
xticks([0 1 2 3 4 5]);
latextitle('\"Ubergangsfunktion des W\"armetauschers');
latexxlabel('$$t$$ in min');
latexylabel('$$T_2$$');
rotateY;


%%    Stochastisch gestörte Übergangsfunktion des Wärmetauschers
%
%
figure(2);
dt=0.025;
Te=3.5;
T=[0:dt:Te];
[Y2, T]=step(Waermetauscher, T);
Y2=10*Y2;     % Temperaturänderung um 10 Kelvin
F=rand(size(Y2))-0.5*ones(size(Y2));
randvalue=0.8;  % Interval der Gleichverteilung [0, randvalue]
Y2r = Y2 + randvalue*F;
Tt=0.5;
kneu=Tt/dt;
Y2t=zeros((Te+Tt)/dt,1);
for k=1:length(T)
    Y2tr(k+kneu,1)=Y2r(k);
    Y2t(k+kneu,1)=Y2(k);
    Tneu(k+kneu)=T(k);
end
Tneu(1:kneu)=[-Tt:dt:-dt];
Y2tr(1:kneu,1)=randvalue*(rand(kneu,1)-0.5*ones(kneu,1));
Y2t(1:kneu,1)=0;
plot(Tneu, Y2tr+57*ones(size(Y2tr)), 'b');  % +57  wegen Arbeitspunktwert 75
axis([-0.5 3.5 56 63]);
yticks([56 57 58 59 60 61 62 63]);
grid on;
%xticks([255 257 259]);
set(gca, 'XTickLabelMode', 'manual');
set(gca, 'XTickLabel', [255:0.5:259]);
latexxlabel('$$t$$ in min');
latexylabel('$$y$$  in $$^o{\rm C}$$');
rotateY;
latextitle('Experimentell ermittelte \"Ubergangsfunktion')


%%              Experimentdaten und Approximation
%
figure(3)
subplot(2,1,1);
plot(Tneu+255, Y2tr+57*ones(size(Y2tr)), 'b', Tneu+255, Y2t+57*ones(size(Y2tr)), 'b',...
    [-5 5]+255, [57 57], 'k:',  [0 0]+255, [0 100], 'k:');
axis([254.5 258 56 63]);
yticks([56 58 60 62]);
xticks([255.0 255.5 256.0 256.5 257.0 257.5 258.0 258.5 259.0]);
latexylabel('$$y$$  in $$^o{\rm C}$$');
rotateY;
latextitle('Analyse der Experimentdaten');
%
subplot(2,1,2);
plot(Ty, Y, 'b', [-0.5 0.8],[0.63*ks, 0.63*ks],'b--',...
    [-0.5 2],[0.95*ks, 0.95*ks],'b--', [-0.5 5], [ks, ks], 'b--',...
    [-0.5 3], [0 0], 'k:');
hold on
plot([-0.5 0],[0 0], 'b');   % Verlängerung der Kurve
plot([0 0], [0 1], 'k:');
latextext(3.6, ks, '$$k_{\rm s}$$');
latextext(0.9, 0.6*ks, '$$0,\!63 k_{\rm s}$$');
latextext(2.2, 0.85*ks, '$$0,\!95 k_{\rm s}$$');
latexylabel('$$y$$ in K');
rotateY;
latexxlabel('$$t$$ in min');
xticks([0  1 2 3]);
axis([-0.5 3 -0.1 0.6])
yticks([0 0.2 0.4 0.6]);
xticks([0 1 2 3 4]);
hold off
%

%%  Vergleich von Simulationsmodell und identifiziertem Modell
%
figure(4)
Modell=ss(-1/0.7, 1/0.7, 0.5, 0);
[Yappr, Tappr]=step(Modell);
%
subplot(4,1,[1,2]);
plot(Ty, Y, 'b--', Tappr, Yappr, 'b',...
      [-1 0],[0 0], 'b');
latexxlabel('$$t$$ in min');
latexylabel('$$y$$ in K');
rotateY;
latextitle('Vergleich von Modell und Experiment');
axis([-0.5 3 -0.1 0.6])
yticks([0 0.2 0.4 0.6]);
xticks([0 1 2 3 4]);


%%     Bilder in .EPS-Dateien
%
Epsbildc15('kennwrt');
