Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
28 commits
Select commit Hold shift + click to select a range
a4c0622
Updated materials streams with Raoult law
biananino May 28, 2024
932568b
Updated Heat exchanger
biananino May 28, 2024
2bf74ee
Updated Mixer
biananino May 28, 2024
0a04041
Updated Separator
biananino May 28, 2024
eabaacf
Updated Shortcut column
biananino May 28, 2024
eb32f47
Updated Flash
biananino May 28, 2024
c9c280c
Updated Splitter
biananino May 28, 2024
1dbc4d3
Updated adabiatic compressor and expander
biananino May 28, 2024
0d7d92e
Updated distillation column components
biananino May 28, 2024
dbdff39
Fixed typo in constraint
biananino May 28, 2024
1d93e3d
Updated main distillation column. Updated example since update is not…
biananino May 28, 2024
456b3c6
Updated conversion reactor
biananino May 28, 2024
e21cf81
Updated NRTL thermodynamic package
biananino May 28, 2024
0ea75cd
Updated material stream NRTL example
biananino May 28, 2024
e10a35e
Updated GraysonStreed and material stream example
biananino May 28, 2024
05588d7
Updated NRTL Binary Phase Envelope
biananino May 28, 2024
47e0c91
Updated PR Binary Phase Envelope
biananino May 28, 2024
f5d5362
Updated UNIFAC Binary Phase Envelope
biananino May 28, 2024
790edbd
Updated UNIQUAC Binary Phase Envelope
biananino May 28, 2024
f20612b
Updated PFR. Doesn't work correctly yet.
biananino May 29, 2024
e4aab4f
Fixed problems with PFR
biananino Jun 14, 2024
892aef7
Updated Absorption Column.
biananino Jun 14, 2024
2ba9e3b
Updated Equilibrium Reactor
biananino Jun 18, 2024
92e085c
Updated UNIQUAC and Material Stream example
biananino Jun 19, 2024
1c037b6
Updated Peng Robinson
biananino Jun 19, 2024
63127c7
Updated UNIFAC
biananino Jun 19, 2024
4b30bdd
Fixed minor mistake in Valve, Centrifugal Pump
biananino Jun 19, 2024
1bef7cd
Fixed minor mistake in Absorbtion Column.
biananino Jun 20, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
70 changes: 38 additions & 32 deletions Simulator/BinaryPhaseEnvelope/NRTL.mo
Original file line number Diff line number Diff line change
Expand Up @@ -4,52 +4,58 @@ package NRTL
extends Modelica.Icons.ExamplesPackage;
model NRTLmodel
import Simulator.Files.ThermodynamicFunctions.*;
gammaNRTLmodel Gamma(NOC = NOC, comp = comp, molFrac = x[:], T = T);
Real density[NOC], BIPS[NOC, NOC, 2];
parameter Integer Nc "Number of components";
parameter Simulator.Files.ChemsepDatabase.GeneralProperties C[Nc] "Component instances array";
gammaNRTLmodel Gamma(Nc = Nc, C = C, molFrac = x[:], T = T);
Real gamma[Nc];
Real x[Nc];
Real T(start = 300);
Real P;
Real K[Nc];
Real density[Nc], BIPS[Nc, Nc, 2];
equation
gamma = Gamma.gamma;
BIPS = Gamma.BIPS;
for i in 1:NOC loop
density[i] = Dens(comp[i].LiqDen, comp[i].Tc, T, P);
for i in 1:Nc loop
density[i] = Dens(C[i].LiqDen, C[i].Tc, T, P);
end for;
for i in 1:NOC loop
K[i] = gamma[i] * Psat(comp[i].VP, T) / P;
for i in 1:Nc loop
K[i] = gamma[i] * Psat(C[i].VP, T) / P;
end for;
end NRTLmodel;

model gammaNRTLmodel
parameter Integer NOC;
parameter Simulator.Files.ChemsepDatabase.GeneralProperties comp[NOC];
Real molFrac[NOC], T;
Real gamma[NOC];
Real tau[NOC, NOC], G[NOC, NOC], alpha[NOC, NOC], A[NOC, NOC], BIPS[NOC, NOC, 2];
Real sum1[NOC], sum2[NOC];
parameter Integer Nc "Number of components";
parameter Simulator.Files.ChemsepDatabase.GeneralProperties C[Nc] "Component instances array";
Real molFrac[Nc], T(start = 300);
Real gamma[Nc];
Real tau[Nc, Nc], G[Nc, Nc], alpha[Nc, Nc], A[Nc, Nc], BIPS[Nc, Nc, 2];
Real sum1[Nc], sum2[Nc];
constant Real R = 1.98721;
equation
A = BIPS[:, :, 1];
alpha = BIPS[:, :, 2];
tau = A ./ (R * T);
for i in 1:NOC loop
for j in 1:NOC loop
for i in 1:Nc loop
for j in 1:Nc loop
G[i, j] = exp(-alpha[i, j] * tau[i, j]);
end for;
end for;
for i in 1:NOC loop
for i in 1:Nc loop
sum1[i] = sum(molFrac[:] .* G[:, i]);
sum2[i] = sum(molFrac[:] .* tau[:, i] .* G[:, i]);
end for;
for i in 1:NOC loop
for i in 1:Nc loop
log(gamma[i]) = sum(molFrac[:] .* tau[:, i] .* G[:, i]) / sum(molFrac[:] .* G[:, i]) + sum(molFrac[:] .* G[i, :] ./ sum1[:] .* (tau[i, :] .- sum2[:] ./ sum1[:]));
end for;
end gammaNRTLmodel;

model base
import data = Simulator.Files.ChemsepDatabase;
parameter Integer NOC;
parameter Real BIP[NOC, NOC, 2];
parameter data.GeneralProperties comp[NOC];
parameter Integer Nc;
parameter Real BIP[Nc, Nc, 2];
parameter Simulator.Files.ChemsepDatabase.GeneralProperties C[Nc] "Component instances array";
extends NRTLmodel(BIPS = BIP);
Real P, T(start = 300), gamma[NOC], K[NOC], x[NOC](each start = 0.5), y[NOC];
Real P, T(start = 300), gamma[Nc], K[Nc], x[Nc](each start = 0.5), y[Nc];
equation
y[:] = K[:] .* x[:];
sum(x[:]) = 1;
Expand All @@ -61,11 +67,11 @@ package NRTL
import data = Simulator.Files.ChemsepDatabase;
parameter data.Onehexene ohex;
parameter data.Acetone acet;
parameter Integer NOC = 2;
parameter Real BIP[NOC, NOC, 2] = Simulator.Files.ThermodynamicFunctions.BIPNRTL(NOC, comp.CAS);
parameter data.GeneralProperties comp[NOC] = {ohex, acet};
base points[41](each P = 1013250, each NOC = NOC, each comp = comp, each BIP = BIP);
Real x[41, NOC], y[41, NOC], T[41];
parameter Integer Nc = 2;
parameter Real BIP[Nc, Nc, 2] = Simulator.Files.ThermodynamicFunctions.BIPNRTL(Nc, C.CAS);
parameter data.GeneralProperties C[Nc] = {ohex, acet};
base points[41](each P = 1013250, each Nc = Nc, each C = C, each BIP = BIP);
Real x[41, Nc], y[41, Nc], T[41];
equation
points[:].x = x;
points[:].y = y;
Expand All @@ -80,11 +86,11 @@ package NRTL
import data = Simulator.Files.ChemsepDatabase;
parameter data.Onehexene ohex;
parameter data.Acetone acet;
parameter Integer NOC = 2;
parameter Real BIP[NOC, NOC, 2] = Simulator.Files.ThermodynamicFunctions.BIPNRTL(NOC, comp.CAS);
parameter data.GeneralProperties comp[NOC] = {ohex, acet};
base points[41](each T = 424, each NOC = NOC, each comp = comp, each BIP = BIP);
Real x[41, NOC], y[41, NOC], P[41];
parameter Integer Nc = 2;
parameter Real BIP[Nc, Nc, 2] = Simulator.Files.ThermodynamicFunctions.BIPNRTL(Nc, C.CAS);
parameter data.GeneralProperties C[Nc] = {ohex, acet};
base points[41](each T = 424, each Nc = Nc, each C = C, each BIP = BIP);
Real x[41, Nc], y[41, Nc], P[41];
equation
points[:].x = x;
points[:].y = y;
Expand All @@ -93,4 +99,4 @@ package NRTL
x[i, 1] = 0 + (i - 1) * 0.025;
end for;
end Pxy;
end NRTL;
end NRTL;
130 changes: 66 additions & 64 deletions Simulator/BinaryPhaseEnvelope/PR.mo
Original file line number Diff line number Diff line change
Expand Up @@ -4,12 +4,12 @@ package PR
extends Modelica.Icons.ExamplesPackage;
function CompresseblityFactor
extends Modelica.Icons.Function;
input Real b[NOC];
input Real aij[NOC, NOC];
input Real b[Nc];
input Real aij[Nc, Nc];
input Real P;
input Real T;
input Integer NOC;
input Real m[NOC];
input Integer Nc;
input Real m[Nc];
output Real am;
output Real bm;
output Real A;
Expand All @@ -20,94 +20,96 @@ package PR
Real C[4];
Real ZR[3, 2];
algorithm
am := sum({{m[i] * m[j] * aij[i, j] for i in 1:NOC} for j in 1:NOC});
am := sum({{m[i] * m[j] * aij[i, j] for i in 1:Nc} for j in 1:Nc});
bm := sum(b .* m);
A := am * P / (R * T) ^ 2;
B := bm * P / (R * T);
C[1] := 1;
C[2] := B - 1;
C[3] := A - 3 * B ^ 2 - 2 * B;
C[4] := B ^ 3 + B ^ 2 - A * B;
ZR := Modelica.Math.Vectors.Utilities.roots(C);
ZR := Modelica.Math.Polynomials.roots(C);
Z := {ZR[i, 1] for i in 1:3};
end CompresseblityFactor;

model PR
parameter Simulator.Files.ChemsepDatabase.GeneralProperties comp[NOC];
parameter Integer NOC;
parameter Simulator.Files.ChemsepDatabase.GeneralProperties C[Nc];
parameter Integer Nc;
parameter Real R = 8.314;
parameter Real kij[NOC, NOC] = Simulator.Files.ThermodynamicFunctions.BIPPR(NOC, comp.name);
Real Tr[NOC];
Real b[NOC];
Real m[NOC];
Real q[NOC];
Real a[NOC];
Real aij[NOC, NOC];
parameter Real kij[Nc, Nc] = Simulator.Files.ThermodynamicFunctions.BIPPR(Nc, C.name);
Real x[Nc];
Real y[Nc];
Real Tr[Nc];
Real b[Nc];
Real m[Nc];
Real q[Nc];
Real a[Nc];
Real aij[Nc, Nc];
Real amL, bmL;
Real AL, BL, Z_L[3];
Real ZL;
Real sum_xa[NOC];
Real liqfugcoeff[NOC];
Real sum_xa[Nc];
Real liqfugcoeff[Nc];
Real amV, bmV;
Real AV, BV, Z_V[3];
Real ZV;
Real sum_ya[NOC];
Real vapfugcoeff[NOC];
Real sum_ya[Nc];
Real vapfugcoeff[Nc];
Real P;
Real T(start = 273);
Real Psat[NOC];
Real Psat[Nc];
//Bubble and Dew Point Calculation
Real Tr_bubl[NOC];
Real a_bubl[NOC];
Real aij_bubl[NOC, NOC];
Real Psat_bubl[NOC];
Real Tr_bubl[Nc];
Real a_bubl[Nc];
Real aij_bubl[Nc, Nc];
Real Psat_bubl[Nc];
Real amL_bubl, bmL_bubl;
Real AL_bubl, BL_bubl, Z_L_bubl[3];
Real ZL_bubl;
Real sum_xa_bubl[NOC];
Real liqfugcoeff_bubl[NOC];
Real gammaBubl[NOC];
Real AL_bubl, BL_bubl(start = 0.005), Z_L_bubl[3];
Real ZL_bubl (each start = 0.005);
Real sum_xa_bubl[Nc];
Real liqfugcoeff_bubl[Nc];
Real gammaBubl[Nc];
Real Tbubl(start = 273);
equation
for i in 1:NOC loop
Psat_bubl[i] = Simulator.Files.ThermodynamicFunctions.Psat(comp[i].VP, Tbubl);
Psat[i] = Simulator.Files.ThermodynamicFunctions.Psat(comp[i].VP, T);
for i in 1:Nc loop
Psat_bubl[i] = Simulator.Files.ThermodynamicFunctions.Psat(C[i].VP, Tbubl);
Psat[i] = Simulator.Files.ThermodynamicFunctions.Psat(C[i].VP, T);
end for;
//Bubble Point and Dew Point Calculation Routine
Tr_bubl = Tbubl ./ comp.Tc;
//Bubble Point and Dew Point Calculation Routine
Tr_bubl = Tbubl ./ C.Tc;
a_bubl = q .* (1 .+ m .* (1 .- sqrt(Tr_bubl))) .^ 2;
aij_bubl = {{(1 - kij[i, j]) * sqrt(a_bubl[i] * a_bubl[j]) for i in 1:NOC} for j in 1:NOC};
(amL_bubl, bmL_bubl, AL_bubl, BL_bubl, Z_L_bubl) = CompresseblityFactor(b, aij_bubl, P, Tbubl, NOC, x[:]);
aij_bubl = {{(1 - kij[i, j]) * sqrt(a_bubl[i] * a_bubl[j]) for i in 1:Nc} for j in 1:Nc};
(amL_bubl, bmL_bubl, AL_bubl, BL_bubl, Z_L_bubl) = CompresseblityFactor(b, aij_bubl, P, Tbubl, Nc, x[:]);
ZL_bubl = min({Z_L_bubl});
sum_xa_bubl = {sum({x[j] * aij_bubl[i, j] for j in 1:NOC}) for i in 1:NOC};
sum_xa_bubl = {sum({x[j] * aij_bubl[i, j] for j in 1:Nc}) for i in 1:Nc};
liqfugcoeff_bubl = exp(AL_bubl / (BL_bubl * sqrt(8)) * log((ZL_bubl + 2.4142135 * BL_bubl) / (ZL_bubl - 0.414213 * BL_bubl)) .* (b / bmL_bubl .- 2 * sum_xa_bubl / amL_bubl) .+ (ZL_bubl - 1) * (b / bmL_bubl) .- log(ZL_bubl - BL_bubl));
liqfugcoeff_bubl[:] = gammaBubl[:] .* P ./ Psat_bubl[:];
P = sum(gammaBubl[:] .* x[:] .* exp(comp[:].VP[2] + comp[:].VP[3] / Tbubl + comp[:].VP[4] * log(Tbubl) + comp[:].VP[5] .* Tbubl .^ comp[:].VP[6]) ./ liqfugcoeff_bubl[:]);
//Calculation of Temperatures at different compositions
Tr = T ./ comp.Tc;
b = 0.0778 * R * comp.Tc ./ comp.Pc;
m = 0.37464 .+ 1.54226 * comp.AF .- 0.26992 * comp.AF .^ 2;
q = 0.45724 * R ^ 2 * comp.Tc .^ 2 ./ comp.Pc;
P = sum(gammaBubl[:] .* x[:] .* exp(C[:].VP[2] + C[:].VP[3] / Tbubl + C[:].VP[4] * log(Tbubl) + C[:].VP[5] .* Tbubl .^ C[:].VP[6]) ./ liqfugcoeff_bubl[:]);
//Calculation of Temperatures at different compositions
Tr = T ./ C.Tc;
b = 0.0778 * R * C.Tc ./ C.Pc;
m = 0.37464 .+ 1.54226 * C.AF .- 0.26992 * C.AF .^ 2;
q = 0.45724 * R ^ 2 *C.Tc.^(2)./ C.Pc;
a = q .* (1 .+ m .* (1 .- sqrt(Tr))) .^ 2;
aij = {{(1 - kij[i, j]) * sqrt(a[i] * a[j]) for i in 1:NOC} for j in 1:NOC};
//Liquid Phase Calculation Routine
(amL, bmL, AL, BL, Z_L) = CompresseblityFactor(b, aij, P, T, NOC, x[:]);
aij = {{(1 - kij[i, j]) * sqrt(a[i] * a[j]) for i in 1:Nc} for j in 1:Nc};
//Liquid Phase Calculation Routine
(amL, bmL, AL, BL, Z_L) = CompresseblityFactor(b, aij, P, T, Nc, x[:]);
ZL = min({Z_L});
sum_xa = {sum({x[j] * aij[i, j] for j in 1:NOC}) for i in 1:NOC};
sum_xa = {sum({x[j] * aij[i, j] for j in 1:Nc}) for i in 1:Nc};
liqfugcoeff = exp(AL / (BL * sqrt(8)) * log((ZL + 2.4142135 * BL) / (ZL - 0.414213 * BL)) .* (b / bmL .- 2 * sum_xa / amL) .+ (ZL - 1) * (b / bmL) .- log(ZL - BL));
//Vapour Phase Calculation Routine
(amV, bmV, AV, BV, Z_V) = CompresseblityFactor(b, aij, P, T, NOC, y[:]);
//Vapour Phase Calculation Routine
(amV, bmV, AV, BV, Z_V) = CompresseblityFactor(b, aij, P, T, Nc, y[:]);
ZV = max({Z_V});
sum_ya = {sum({y[j] * aij[i, j] for j in 1:NOC}) for i in 1:NOC};
sum_ya = {sum({y[j] * aij[i, j] for j in 1:Nc}) for i in 1:Nc};
vapfugcoeff = exp(AV / (BV * sqrt(8)) * log((ZV + 2.4142135 * BV) / (ZV - 0.414213 * BV)) .* (b / bmV .- 2 * sum_ya / amV) .+ (ZV - 1) * (b / bmV) .- log(ZV - BV));
end PR;

model PhaseEquilibria
import data = Simulator.Files.ChemsepDatabase;
parameter data.Ethane eth;
parameter data.Propane prop;
extends PR(NOC = 2, comp = {eth, prop});
Real P, T(start = 273), K[NOC], x[NOC](each start = 0.5), y[NOC], Tbubl(start = 273);
extends PR(Nc = 2, C = {eth, prop});
Real P, T(start = 273), K[Nc], x[Nc](each start = 0.5), y[Nc], Tbubl(start = 273);
equation
K[:] = liqfugcoeff[:] ./ vapfugcoeff[:];
y[:] = K[:] .* x[:];
Expand All @@ -120,13 +122,13 @@ package PR
import data = Simulator.Files.ChemsepDatabase;
parameter data.Ethane eth;
parameter data.Propane prop;
parameter Integer NOC = 2;
parameter Integer N = 2;
parameter data.GeneralProperties comp[NOC] = {eth, prop};
PhaseEquilibria points[N](each T = 210, each NOC = NOC, each comp = comp, each T(start = 273), each Tbubl(start = 273), each x(each start = 0.5), each y(each start = 0.5));
parameter Integer Nc = 2;
parameter Integer N = 10;
parameter data.GeneralProperties C[Nc] = {eth, prop};
PhaseEquilibria points[N](each T = 210, each Nc = Nc, each C = C, each T(start = 273), each Tbubl(start = 273), each x(each start = 0.5), each y(each start = 0.5));
Real x1[N], y1[N], x2[N], y2[N], P[N](each start = 101325), Tbubl[N], Temp[N];
equation
//Generation of Points to compute Bubble Temperature
//Generation of Points to Compute Bubble Temperature
points[:].x[1] = x1[:];
points[:].y[1] = y1[:];
points[:].x[2] = x2[:];
Expand All @@ -148,11 +150,11 @@ package PR
import data = Simulator.Files.ChemsepDatabase;
parameter data.Ethane eth;
parameter data.Propane prop;
parameter Integer NOC = 2;
parameter Integer N = 10;
parameter data.GeneralProperties comp[NOC] = {eth, prop};
PhaseEquilibria points[N](each P = 101325, each NOC = NOC, each comp = comp, each T(start = 273), each Tbubl(start = 273), each x(each start = 0.5), each y(each start = 0.5));
Real x[N, NOC], y[N, NOC], T[N], Tbubl[N], T_PR[N];
parameter Integer Nc = 2;
parameter Integer N = 21;
parameter data.GeneralProperties C[Nc] = {eth, prop};
PhaseEquilibria points[N](each P = 101325, each Nc = Nc, each C = C, each T(start = 273), each Tbubl(start = 273), each x(each start = 0.5), each y(each start = 0.5));
Real x[N, Nc], y[N, Nc], T[N], Tbubl[N], T_PR[N];
equation
points[:].x = x;
points[:].y = y;
Expand All @@ -164,7 +166,7 @@ package PR
T_PR[i] = T[i];
end for;
for i in 1:N loop
x[i, 1] = 0 + (i - 1) * 0.025;
x[i, 1] = 0.5 + (i - 1) * 0.025;
end for;
end Txy;
end PR;
end PR;
8 changes: 4 additions & 4 deletions Simulator/BinaryPhaseEnvelope/UNIFAC.mo
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ package UNIFAC
//Libraries
import Simulator.*;
//Extension of Chemsep Database
Simulator.Files.ChemsepDatabase data;
import data = Simulator.Files.ChemsepDatabase;
//Parameter Section
//Selection of compounds
parameter data.Methylethylketone meth;
Expand Down Expand Up @@ -135,7 +135,7 @@ package UNIFAC
//Libraries
import Simulator.*;
//Extension of Chemsep Database
Simulator.Files.ChemsepDatabase data;
import data = Simulator.Files.ChemsepDatabase;
//Parameter Section
//Selection of compounds
parameter data.Methylethylketone meth;
Expand Down Expand Up @@ -177,7 +177,7 @@ package UNIFAC
//Mole Fractions (x-axis) of the T-x-y plot
Real z1[N + 1], z2[N + 1];
//Bubble Temperature
Real T[N + 1](unit = "K", each start = 300);
Real T[N + 1](each unit = "K", each start = 300);
//Distribution coefficient
Real K1[N + 1];
//Vapour Phase Mole Fraction
Expand Down Expand Up @@ -264,4 +264,4 @@ package UNIFAC
end Txy;

//================================================================================================================
end UNIFAC;
end UNIFAC;
Loading