REQUIRE "thermodynamics.a4l";

ATOM virialC REFINES solver_var
DIMENSION L^6/Q^2
DEFAULT 7000.0 {cm^6/mol^2};
lower_bound := 0.0 {cm^6/mol^2};
upper_bound := 1.0 {m^6/mol^2};
nominal := 7000.0 {cm^6/mol^2};
END virialC;

ATOM virialB REFINES solver_var
DIMENSION L^3/Q
DEFAULT 1.0 {m^3/mol};
lower_bound := -1000.00 {m^3/mol};
upper_bound := 1000.0 {m^3/mol};
nominal := 1.0 {m^3/mol};
END virialB;

MODEL parta REFINES cmumodel;

R IS_A gas_constant;
nc IS_A set OF symbol_constant;
nc :== ['ethylene'];
cd IS_A components_data(['ethylene'], 'ethylene');
P, Pc IS_A pressure;
T, Tc IS_A temperature;
T_degC, Z IS_A factor;
pvc IS_A Pitzer_vapor_component(P, T, cd.data['ethylene']);
Tr, Pr IS_A factor;
B, V IS_A molar_volume;
C IS_A virialC;

Pc = pvc.data.Pc;
Tc = pvc.data.Tc;

eq1: Tr = T/Tc;
eq2: Pr = P/Pc;
eq3: T_degC = T/1{K} - 273.15;

(* for part a *)
eq4: Z*R*T = P*V;
eq5: Z*V^2 = V^2 + B*V + C;
(* for part a *)


METHODS
        METHOD doall;
        RUN ClearAll;
        P.fixed := TRUE;
        P := 12.0 {bar};
        T_degC.fixed := TRUE;
        T_degC := 25.0;
        B.fixed := TRUE;
        B := -140.0 {cm^3/mol};
        C.fixed := TRUE;
        C := 7200.0 {cm^6/mol^2};
        END doall;

END parta;

MODEL partb REFINES cmumodel;

R IS_A gas_constant;
nc IS_A set OF symbol_constant;
nc :== ['ethylene'];
cd IS_A components_data(['ethylene'], 'ethylene');
P, Pc IS_A pressure;
T, Tc IS_A temperature;
T_degC, Z IS_A factor;
pvc IS_A Pitzer_vapor_component(P, T, cd.data['ethylene']);
Tr, Pr, omega IS_A factor;
V IS_A molar_volume;
B IS_A virialB;
Bhat, B0, B1 IS_A factor;

Pc = pvc.data.Pc;
Tc = pvc.data.Tc;
omega = pvc.data.omega;

eq1: Tr = T/Tc;
eq2: Pr = P/Pc;
eq3: T_degC = T/1{K} - 273.15;

eq4: Bhat = B0 + omega*B1;
eq5: B0*Tr^(1.6) = 0.083*Tr^(1.6) - 0.422;
eq6: B1*Tr^(4.2) = 0.139*Tr^(4.2) - 0.172;

eq7: Z = 1.0 + B0*Pr/Tr + omega*B1*Pr/Tr;

eq8: Z*R*T = P*V;
(*eq9: Z = 1.0 + B*P/(R*T);   *)
eq10: B = Bhat*1.0{m^3/mol};


METHODS
        METHOD doall;
        RUN ClearAll;
        P.fixed := TRUE;
        P := 12.0 {bar};
        T_degC.fixed := TRUE;
        T_degC := 25.0;
        END doall;

END partb;

MODEL partc REFINES cmumodel;
(* rk *)
R IS_A gas_constant;
nc IS_A set OF symbol_constant;
nc :== ['ethylene'];
cd IS_A components_data(['ethylene'], 'ethylene');
P, Pc IS_A pressure;
T, Tc IS_A temperature;
alpha, T_degC, Z IS_A factor;
pvc IS_A Pitzer_vapor_component(P, T, cd.data['ethylene']);
eps, sigma, beta, q,Tr, Pr, omega, OMEGA, PSI IS_A factor;
V IS_A molar_volume;

Pc = pvc.data.Pc;
Tc = pvc.data.Tc;
omega = pvc.data.omega;

eq1: Tr = T/Tc;
eq2: Pr = P/Pc;
eq3: T_degC = T/1{K} - 273.15;

eq4: alpha = Tr^(0.5);
eq5: q = PSI*alpha/(OMEGA*Tr);
eq6: beta = OMEGA*Pr/Tr;
eq7: Z = 1.0 + beta - q*beta*(Z - beta)/((Z + eps*beta)*(Z + sigma*beta));
eq8: P*V = Z*R*T;
METHODS
        METHOD doall;
        RUN ClearAll;
        P.fixed := TRUE;
        P := 12.0 {bar};
        T_degC.fixed := TRUE;
        T_degC := 25.0;
        sigma.fixed := TRUE;
        sigma := 1.0;
        eps.fixed := TRUE;
        eps := 0.0;
        OMEGA.fixed := TRUE;
        OMEGA := 0.08664;
        PSI.fixed := TRUE;
        PSI := 0.42748;
        Z.lower_bound := 0.0;
        Z.upper_bound := 10.0;
        Z.nominal := 1.0;
        END doall;

END partc;

MODEL partd REFINES cmumodel;
(* srk *)
R IS_A gas_constant;
nc IS_A set OF symbol_constant;
nc :== ['ethylene'];
cd IS_A components_data(['ethylene'], 'ethylene');
P, Pc IS_A pressure;
T, Tc IS_A temperature;
alpha, T_degC, Z IS_A factor;
pvc IS_A Pitzer_vapor_component(P, T, cd.data['ethylene']);
eps, sigma, beta, q,Tr, Pr, omega, OMEGA, PSI IS_A factor;
V IS_A molar_volume;

Pc = pvc.data.Pc;
Tc = pvc.data.Tc;
omega = pvc.data.omega;

eq1: Tr = T/Tc;
eq2: Pr = P/Pc;
eq3: T_degC = T/1{K} - 273.15;

eq4: alpha = (1.0 + (0.480+1.574*omega-0.176*omega^2)*(1.0-Tr^(0.5)))^2;
eq5: q = PSI*alpha/(OMEGA*Tr);
eq6: beta = OMEGA*Pr/Tr;
eq7: Z = 1.0 + beta - q*beta*(Z - beta)/((Z + eps*beta)*(Z + sigma*beta));
eq8: P*V = Z*R*T;
METHODS
        METHOD doall;
        RUN ClearAll;
        P.fixed := TRUE;
        P := 12.0 {bar};
        T_degC.fixed := TRUE;
        T_degC := 25.0;
        sigma.fixed := TRUE;
        sigma := 1.0;
        eps.fixed := TRUE;
        eps := 0.0;
        OMEGA.fixed := TRUE;
        OMEGA := 0.0864;
        PSI.fixed := TRUE;
        PSI := 0.42748;
        Z.lower_bound := 0.0;
        Z.upper_bound := 10.0;
        Z.nominal := 1.0;
        END doall;

END partd;

MODEL parte REFINES cmumodel;

R IS_A gas_constant;
nc IS_A set OF symbol_constant;
nc :== ['ethylene'];
cd IS_A components_data(['ethylene'], 'ethylene');
P, Pc IS_A pressure;
T, Tc IS_A temperature;
alpha, T_degC, Z IS_A factor;
pvc IS_A Pitzer_vapor_component(P, T, cd.data['ethylene']);
eps, sigma, beta, q,Tr, Pr, omega, OMEGA, PSI IS_A factor;
V IS_A molar_volume;

Pc = pvc.data.Pc;
Tc = pvc.data.Tc;
omega = pvc.data.omega;

eq1: Tr = T/Tc;
eq2: Pr = P/Pc;
eq3: T_degC = T/1{K} - 273.15;

eq4: alpha = (1.0 + (0.37464+1.54226*omega-0.26992*omega^2)*(1.0-Tr^(0.5)))^2;
eq5: q = PSI*alpha/(OMEGA*Tr);
eq6: beta = OMEGA*Pr/Tr;
eq7: Z = 1.0 + beta - q*beta*(Z - beta)/((Z + eps*beta)*(Z + sigma*beta));
eq8: P*V = Z*R*T;
METHODS
        METHOD doall;
        RUN ClearAll;
        P.fixed := TRUE;
        P := 12.0 {bar};
        T_degC.fixed := TRUE;
        T_degC := 25.0;
        sigma.fixed := TRUE;
        sigma := 2.41424;
        eps.fixed := TRUE;
        eps := -0.41424;
        OMEGA.fixed := TRUE;
        OMEGA := 0.07779;
        PSI.fixed := TRUE;
        PSI := 0.45724;
        Z.lower_bound := -50.0;
        Z.upper_bound := 100.0;
        Z.nominal := 0.05;
        END doall;

END parte;