% FordamperOver.m
% Finner entalpidifferansen over fordamper
% Utgangspunkt: 3P CO2 SC IC
% h = T*s + my*N
% 
% Startverdier som funker: T1=8,-10, --> V1=5e-5
%                          T1=-20    --> V1=4e-5
clear all

p1 = 48.30e5;           % [Pa] Gitt
T1 = 283;         % [K]   "
T2 = T1-21;             % [K]  Konst T ved fordampning
overheting = 0;
V1 = 4e-5;              % [m3] Startverdi, væskevolum
N = [1];                % En komponent (ren co2)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Overheting

% PENG-ROBINSON
% Før fordamper:
v1  = V1;
dv1 = v1;
% Newtoniterasjoner til V ved ønsket trykk finnes:
while abs(dv1)>1e-6*v1
    a1 = co2_pr([T1,v1,N]');
    dv1 = (p1 + a1.g(2))/(-a1.H(2,2));
    v1 = v1 + dv1;
end
disp('PR TILFELLE B: Overhetet damp')
h1 =-a1.g(1)*T1 + a1.g(3)*N;   % Entalpi før fordamper
disp(['Entalpi før fordamper:   H1 = ',num2str(h1),' J/mol'])

% Etter fordamper, tørrmettet:
pH_co2_pr;
v2  = V;
a2  = co2_pr([T2,v2,N]');
p2 = -a2.g(2);
% Overheting, konst p2:
t2 = T2 + overheting;      % [K] Overhetet damp (Tevap=ca -10C ved 25bar)
dv2= v2;
while abs(dv2)>1e-6*v2
    a2 = co2_pr([t2,v2,N]');
    dv2 = (p2 + a2.g(2))/(-a2.H(2,2));
    v2 = v2 + dv2;
end
h2 = -a2.g(1)*t2 + a2.g(3)*N;
disp(['Entalpi etter fordamper: H2 = ',num2str(h2),' J/mol'])
% Differanse:
Dh = h2 - h1;
disp(['Differanse:           H2-H1 =   ',num2str(Dh),' J/mol'])
disp(' ')


% SPAN-WAGNER
% Før fordamper:
dV1 = V1;
while abs(dV1)>1e-6*V1
    A1 = co2_sw([T1,V1,N]');
    dV1 = (p1 + A1.g(2))/(-A1.H(2,2));
    V1 = V1 + dV1;
end
disp('SW TILFELLE B: Overhetet damp')
H1 =-A1.g(1)*T1 + A1.g(3)*N;   % Entalpi før fordamper
disp(['Entalpi før fordamper:   H1 = ',num2str(H1),' J/mol'])

% Etter fordamper, tørrmettet:
pH_co2_sw;
V2 = V;
A2  = co2_sw([T2,V2,N]');
P2 = -A2.g(2);
% Overheting:
T2 = T2 + overheting;
dV2= V2;
while abs(dV2)>1e-6*V2
    A2 = co2_sw([T2,V2,N]');
    dV2 = (P2 + A2.g(2))/(-A2.H(2,2));
    V2 = V2 + dV2;
end
H2 = -A2.g(1)*T2 + A2.g(3)*N;
disp(['Entalpi etter fordamper: H2 = ',num2str(H2),' J/mol'])

% Differanse:
DH = H2 - H1;
disp(['Differanse:           H2-H1 =   ',num2str(DH),' J/mol'])
disp(' ')


