GOAT/calcab_multicomp.m at main · smithh2710/GOAT · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
%% CALCULATE DIMENSIONLESS ATTRACTION AND COVOLUME, A & B
% -------------------------------------------------------------------------
% INPUTS:
%   press      - pressure [Pa]
%   temp       - temperature [K]
%   pressc     - critical pressure [Pa]
%   tempc      - critical temperature [K]
%   acentric   - acentric factor
%   eos_type   - 'PR' (default) or 'SRK'
%   components - (optional) cell array of component names for tc-SRK lookup
%
% OUTPUTS:
%   A, B       - Dimensionless EOS parameters
%
% For SRK: Uses Pina-Martinez (2018) Twu91 alpha with component-specific
%          L, M, N from get_tc_srk_parameters(). Falls back to generalized
%          correlation if component not in table.

function [A, B] = calcab_multicomp(press, temp, pressc, tempc, acentric, eos_type, components)

if nargin < 6 || isempty(eos_type)
    eos_type = 'PR';
end

if nargin < 7
    components = {};
end

ncomp = size(pressc,1);
m = zeros(ncomp,1);
alpha = zeros(ncomp,1);
A = zeros(ncomp,1);
B = zeros(ncomp,1);
L = zeros(ncomp,1);
M = zeros(ncomp,1);
N = zeros(ncomp,1);

switch upper(eos_type)
    case 'PR'
        omegaa = 0.45724;
        omegab = 0.0778;
    case 'SRK'
        omegaa = 0.427480012178;
        omegab = 0.08664000034332;
end

for i = 1:ncomp

    pressr = press/pressc(i);
    tempr = temp/tempc(i);

    switch upper(eos_type)
        case 'PR'
            if acentric(i) > 0.49
                m(i) = 0.379642 + 1.48503*acentric(i) - 0.164423*acentric(i)^2 + 0.016666*acentric(i)^3;
            else
                m(i) = 0.37464 + 1.54226*acentric(i) - 0.26992*acentric(i)^2;
             end
             alpha(i) = ( 1 + m(i)*(1 - sqrt(tempr)) )^2;

            % m(i) = 0.384401 + 1.52276 * acentric(i)-0.213808* acentric(i)^2 + 0.034616 * acentric(i)^3 - 0.001976 * acentric(i)^4 ; % m as suggested by magoulus and tassios

            % pina-martinez(2019) alpha rule
            % alpha(i) = (1 + (0.3919 + 1.4996 * acentric(i) - 0.2721 * acentric(i)^2 + 0.1063*acentric(i)^3 ) * (1-tempr^(0.5))  )^2;


            % pina-martinez(2018) alpha rule
            % L(i) = 0.0728 + 0.6693 * acentric(i) + 0.0925 * acentric(i)^2 ;
            % M(i) = 0.8788 - 0.2258* acentric(i) + 0.1695* acentric(i)^2 ;
            % alpha(i) = tempr^(2* (M(i)-1)) * exp(L(i)*(1-tempr^(2*M(i)))) ;
            %
            % alpha li-yang (2012)
            % alpha(i)= exp ((0.13280-0.05052* acentric(i) + 0.25948* acentric(i)^2)* (1-tempr) + 0.81769 * log(    1+(0.31355 + 1.86745 *acentric(i)- 0.52604 * acentric(i)^2) * (1-tempr^(0.5))       )^2         ) ;


            % Gasem et al. (2001) alpha function for PR EOS

            % Aa = 2.0; Ba = 0.836; Ca = 0.134; Da = 0.508; Ea = -0.0467;
            % alpha(i) = exp((Aa + Ba*tempr) * (1 - tempr^(Ca + Da*acentric(i) + Ea*acentric(i)^2)));


        case 'SRK'
            found = false;

            % if i <= length(components) && ~isempty(components{i})
            %     [params, found] = get_tc_srk_parameters(components{i});
            %     if found
            %         L(i) = params.L;
            %         M(i) = params.M;
            %         N(i) = params.N;
            %     end
            % end

            if ~found
                L(i) = 0.1359 + 0.7535*acentric(i) + 0.0611*acentric(i)^2;
                M(i) = 0.8787 - 0.2063*acentric(i) + 0.1709*acentric(i)^2;
                N(i) = 2.0;
            end

             alpha(i) = tempr^(N(i)*(M(i)-1)) * exp(L(i)*(1-tempr^(N(i)*M(i))));

            % alpha function by Soave (1972)
              % m(i) = 0.480 + 1.574*acentric(i) - 0.176*acentric(i)^2;

            % Graboski & Daubert (1978) - improved for heavier hydrocarbons
              % m(i) = 0.48508 + 1.55171*acentric(i) - 0.15613*acentric(i)^2;
              %
                % alpha(i) = ( 1 + m(i)*(1 - sqrt(tempr)) )^2;
    end

    A(i) = omegaa*alpha(i)*pressr/tempr^2;
    B(i) = omegab*pressr/tempr;

end

end