diff --git a/Simulator/BinaryPhaseEnvelope/NRTL.mo b/Simulator/BinaryPhaseEnvelope/NRTL.mo index cc33480..80622a3 100644 --- a/Simulator/BinaryPhaseEnvelope/NRTL.mo +++ b/Simulator/BinaryPhaseEnvelope/NRTL.mo @@ -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; @@ -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; @@ -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; @@ -93,4 +99,4 @@ package NRTL x[i, 1] = 0 + (i - 1) * 0.025; end for; end Pxy; -end NRTL; +end NRTL; \ No newline at end of file diff --git a/Simulator/BinaryPhaseEnvelope/PR.mo b/Simulator/BinaryPhaseEnvelope/PR.mo index f9fa542..ac646fa 100644 --- a/Simulator/BinaryPhaseEnvelope/PR.mo +++ b/Simulator/BinaryPhaseEnvelope/PR.mo @@ -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; @@ -20,7 +20,7 @@ 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); @@ -28,77 +28,79 @@ package PR 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; @@ -106,8 +108,8 @@ package PR 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[:]; @@ -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[:]; @@ -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; @@ -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; \ No newline at end of file diff --git a/Simulator/BinaryPhaseEnvelope/UNIFAC.mo b/Simulator/BinaryPhaseEnvelope/UNIFAC.mo index 59d3bc5..cd767ac 100644 --- a/Simulator/BinaryPhaseEnvelope/UNIFAC.mo +++ b/Simulator/BinaryPhaseEnvelope/UNIFAC.mo @@ -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; @@ -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; @@ -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 @@ -264,4 +264,4 @@ package UNIFAC end Txy; //================================================================================================================ -end UNIFAC; +end UNIFAC; \ No newline at end of file diff --git a/Simulator/BinaryPhaseEnvelope/UNIQUAC.mo b/Simulator/BinaryPhaseEnvelope/UNIQUAC.mo index bb9367f..600563d 100644 --- a/Simulator/BinaryPhaseEnvelope/UNIQUAC.mo +++ b/Simulator/BinaryPhaseEnvelope/UNIQUAC.mo @@ -14,12 +14,12 @@ package UNIQUAC input Real R[NOC], Q[NOC]; input Real tow[NOC, NOC]; input Real towk[N + 1, NOC, NOC]; - parameter Real Z = 10 "Compresseblity Factor"; - parameter Real R_gas = 1.98721 "Gas Constant"; //Activity coefficients output Real gammaBubl1[N + 1], gammaBubl2[N + 1]; protected //Intermediate parameters used to calculate the Combinatorial and Residual contribution" + parameter Real Z = 10 "Compresseblity Factor"; + parameter Real R_gas = 1.98721 "Gas Constant"; Real r_bubl[N + 1], q_bubl[N + 1]; Real teta1_bubl[N + 1], teta2_bubl[N + 1]; Real S1_bubl[N + 1], S2_bubl[N + 1]; @@ -87,7 +87,7 @@ package UNIQUAC //Libraries import Simulator.*; //Extension of Chemsep Database - Simulator.Files.ChemsepDatabase data; + import data = Simulator.Files.ChemsepDatabase; //Parameter Section //Selection of compounds parameter data.Water wat; @@ -120,7 +120,7 @@ package UNIQUAC //Vapour Phase Mole Fraction Real y1[N + 1](each start = 0.5), y2[N + 1](each start = 0.5); //Vapour Pressure at the chosen temperature - Real Psat[NOC](unit = "Pa") "Vapour Pressure"; + Real Psat[NOC](each unit = "Pa") "Vapour Pressure"; //========================================================================================= //Equation Section equation @@ -170,7 +170,7 @@ package UNIQUAC //Libraries import Simulator.*; //Extension of Chemsep database - Simulator.Files.ChemsepDatabase data; + import data = Simulator.Files.ChemsepDatabase; //Parameter Section //Selection of compounds parameter data.Water wat; @@ -255,4 +255,4 @@ package UNIQUAC //================================================================================================ //============================================================================================================== //================================================================================================================ -end UNIQUAC; +end UNIQUAC; \ No newline at end of file diff --git a/Simulator/Examples/Absorption.mo b/Simulator/Examples/Absorption.mo index 901a9c8..d94419e 100644 --- a/Simulator/Examples/Absorption.mo +++ b/Simulator/Examples/Absorption.mo @@ -17,8 +17,8 @@ package Absorption "Example of Simulating an Absorption Column" end Tray; model AbsColumn "Extension of Absorption COlumn with instance of Tray" - extends Simulator.UnitOperations.AbsorptionColumn.AbsCol; - Tray tray[Nt](each Nc = Nc, each C = C); + extends Simulator.UnitOperations.AbsorptionColumn.AbsCol( + redeclare Tray tray[Nt](each Nc = Nc, each C = C)); annotation( Documentation(info = "
This is a non-executable model is created inside the package Absorption to extend the Absorption Column model along with the necessary property method from ThermodynamicPackages which is RaoultsLaw in this case.