Commit 27fb292e authored by Andreas Michael Bremen's avatar Andreas Michael Bremen
Browse files

- cleanup of SolidLiquidPhase with some variable name changes for...

- cleanup of SolidLiquidPhase with some variable name changes for initialization (stoichiometry matrices, nullspace etc.)
parent d8109f25
within ElectrolyteMedia.Media.Common.Reaction;
function calc_lambda_id
"calculates lambda with nX*nX identity matrix"
input Real[:,:] lambda;
output Real[size(lambda,1),size(lambda,2)] lambda_id;
output Integer[size(lambda,1),size(lambda,1)] P_id;
protected
Integer nF = size(lambda,1);
Integer nX = size(lambda,2);
Real[nF,nX] LU;
Integer[nX] pivots_LAPACK;
Real[nF,nX] L;
Integer[nF,nF] P;
Real[nX,nX] Lsquare;
Real[nX,nX] LsquareT;
Real[nX,nX] LU_;
Integer[nX] pivots_LAPACK_;
Real[nX,nX] U_;
input Real[:,:] nu_id;
output Real[size(nu_id,2),size(nu_id,2)-size(nu_id,1)] lambda_id;
algorithm
(LU,pivots_LAPACK) :=Modelica.Math.Matrices.LU(lambda);
L :=Media.Common.Reaction.calc_L(LU);
P :=Media.Common.Reaction.calc_P(pivots_LAPACK,nF);
Lsquare :=L[1:nX, :];
LsquareT :=transpose(Lsquare);
(LU_,pivots_LAPACK_) :=Modelica.Math.Matrices.LU(LsquareT);
U_ :=Media.Common.Reaction.calc_U(LU_);
lambda_id :=L*Modelica.Math.Matrices.inv(transpose(U_));
P_id :=P;
lambda_id :=transpose(cat(2,-transpose(nu_id[:, size(nu_id,1) + 1:size(nu_id,2)]),identity(size(nu_id,2)-size(nu_id,1))));
end calc_lambda_id;
......@@ -5,4 +5,3 @@ calc_L
calc_U
calc_P_nu_id
calc_lambda_id
calc_P_lambda_id
within ElectrolyteMedia.Media.LiquidPhase.Common;
record LiquidInteractionDataRecord "Coefficient data record for interaction properties between aqueous solute species"
parameter Integer nLi;
parameter Integer nLi = 1;
//Bromley parameter
Real[nLi,nLi] Bromley_ij = zeros(nLi,nLi) "Bromley parameter";
......
......@@ -613,10 +613,6 @@ package Common
Temperature(min=273.15, max=573.15, start=300, nominal=300),
T_default = userInterface.Tstart,
p_default = userInterface.pstart);
// nXi = nX-1,
// fixedX=false,
// reducedX = true,
// singleState=false,
redeclare record extends ThermodynamicState "Thermodynamic state variables"
Real[nF] Xfull;
......@@ -646,7 +642,7 @@ package Common
//Reaction
parameter Integer nR = 1 "Number of solid-liquid and dissociation equilibria" annotation(Dialog(tab="Reaction"));
Real[nR,ns+nL] nu = zeros(nR,ns+nL) "Stoichiometry matrix of solid-liquid and dissociation equilibria" annotation(Dialog(tab="Reaction"));
parameter Real[nR,ns+nL] nu = cat(2,identity(nR),-ones(nR,ns+nL-nR)) "Stoichiometry matrix of solid-liquid and dissociation equilibria" annotation(Dialog(tab="Reaction"));
//General
String mediumName = "mediumName" "Name of medium" annotation(Dialog(tab="General"));
......@@ -723,35 +719,19 @@ package Common
constant Media.Common.Types.LiquidModel LiquidModel=userInterface.LiquidModel;
constant Real[nR,nF] nu = userInterface.nu;
constant Real[nR,nF] nu_id_= Media.Common.Reaction.calc_nu_id(nu);
//nu
constant Integer[nF,nF] P_to_orig_=Media.Common.Reaction.calc_P_nu_id(nu);
//identity(nF);
constant Integer[nF,nF] P_to_id_ = transpose(P_to_orig_);
constant Real[nR,nF] nu_id__ = nu_id_*transpose(P_to_orig_) "nu in original order";
constant Real[nF,nX] lambda_nu = transpose(cat(2,-transpose(nu_id_[:,nR+1:nF]),identity(nX)));
// constant Real[nF,nX] lambda_id = cat(2,identity(nX),transpose(nu_id[1:nX,nX+1:nF]));
constant Real[nR,nF] nu_mass = Functions.Reaction.calc_nu_mass(nu);
constant Real[nF] MMX_ = P_to_id_*MMX;
constant Real[nF,nX] lambda= Media.Common.Reaction.calc_lambda(nu, nR, nF) "Nullspace of stoichiometry matrix";
constant Real[nF,nX] lambda_mass={{lambda[i,j]/MMX[i] for j in 1:nX} for i in 1:nF} "Nullspace of mass based stoichiometry matrix";
constant Real[nR,nF] nu_id= Media.Common.Reaction.calc_nu_id(nu) "Identity transformation of stoichiometry matrix for initialization";
constant Integer[nF,nF] P_to_orig=Media.Common.Reaction.calc_P_nu_id(nu) "Permutation matrix between original order and identity transformation";
constant Integer[nF,nF] P_to_id = transpose(P_to_orig) "Permutation matrix between identity transformation and original order";
constant Real[nF,nX] lambda_id = Media.Common.Reaction.calc_lambda_id(nu_id) "lambda in identity transformation";// transpose(cat(2,-transpose(nu_id[:,nR+1:nF]),identity(nX)))
constant Real[nF] MMX_id = P_to_id*MMX;
constant UserInterface
userInterface "User interface"
annotation (Placement(transformation(extent={{-10,-8},{10,12}})));
constant Real[nF,nX] lambda_=
Media.Common.Reaction.calc_lambda(
nu,
nR,
nF);
constant Real[nF,nX] lambda_id = Media.Common.Reaction.calc_lambda_id(lambda_);
constant Real[nF,nX] lambda = P_to_orig*lambda_id;
constant Integer[nF,nF] P_to_orig= Media.Common.Reaction.calc_P_lambda_id(lambda_);
constant Integer[nF,nF] P_to_id = transpose(P_to_orig);
constant Real[nR,nF] nu_id = nu*transpose(P_to_id);
constant Real[nF,nX] lambda_mass_=P_to_orig_*{{lambda_nu[i,j]/MMX_[i] for j in 1:nX} for i in 1:nF} "Null space of mass based stoichiometry matrix";
constant Real[nF,nX] lambda_mass={{lambda_[i,j]/MMX[i] for j in 1:nX} for i in 1:nF} "Null space of mass based stoichiometry matrix";
// constant Real[nX] scale = {sum(lambda_mass_[i,j] for i in 1:nF) for j in 1:nX};
// constant Real[:,:] lambda_mass = {{lambda_mass_[i,j]/scale[j] for j in 1:nX} for i in 1:nF};
constant MolarMass MH2O = IF97.MH2O "Molar mass of solvent";
constant SpecificHeatCapacity RH2O = IF97.RH2O;
......@@ -805,7 +785,7 @@ package Common
Tstart,
Xredstart);
parameter MassFraction[nX] Xredstart = reference_X;
parameter Real[nX] Xredstart_ = transpose(lambda_mass_)*Xfullstart "for scaling of species moles consistent with maxx invariant vector X";
parameter Real[nX] Xredstart_ = transpose(lambda_mass)*Xfullstart "for scaling of species moles consistent with maxx invariant vector X";
parameter Real[ns+nL] mfullstart = Xfullstart/sum(Xredstart_);
final parameter Real[nL-1] mstart_ = Functions.LiquidFunctions.calc_mfromX(Xlstart);
final parameter Real[nL] mstart = cat(1,mstart_,{Ylstart[nL]/MH2O});
......@@ -1231,7 +1211,7 @@ package Common
Integer index;
Real[nX] Xs;
Real[nF,nX] lambda_mass_id = P_to_id_ * lambda_mass;
Real[nF,nX] lambda_mass_id = P_to_id * lambda_mass;
algorithm
//ensure no species amount to be zero for numerical issues
for i in 1:nX loop
......@@ -1241,14 +1221,14 @@ package Common
//convert reduced mass fraction basis
Xs :=Modelica.Math.Matrices.solve(transpose(lambda_mass_id[nR+1:nF, :]), Xrednonzero);
Xfull :=cat(1,zeros(nF - nX),Xs);
Xfull :=P_to_orig_*Xfull;
Xfull :=P_to_orig*Xfull;
Xred_ :=transpose(lambda_mass_)*Xfull;
Xred_ :=transpose(lambda_mass)*Xfull;
// find initial guess x1 for Newton solver
if sum(Xinit)> 0 then
// ninit from initial guess provided as input (Xinit)
scale :=sum(transpose(lambda_mass_)*Xinit);
scale :=sum(transpose(lambda_mass)*Xinit);
minit :=Xinit*scale;
ninit :=minit./MMX;
for i in 1:nF loop
......@@ -1263,7 +1243,7 @@ package Common
// Newton algorithm
solutionfound :=false;
x1 :=P_to_id_*x1;
x1 :=P_to_id*x1;
while not solutionfound loop
f :=Initialization.calc_f.SLE_Tp(x1,Xred_,T,p);
J :=Initialization.calc_J.SLE_Tp(x1,T,p);
......@@ -1275,7 +1255,7 @@ package Common
end if;
x1 :=x2;
end while;
x2 :=P_to_orig_*x2;
x2 :=P_to_orig*x2;
Yfull :=x2[1:nF]/sum(x2[1:nF]);
Xfull :=moleToMassFractions(Yfull, MMX);
......@@ -1511,7 +1491,7 @@ package Common
//assert(sum(x[1:ns+nL]) > 0, "x is zero in SLE_Tp");
x_orig :=P_to_orig_*x;
x_orig :=P_to_orig*x;
for i in 1:ns+nL loop
x_orig[i] :=x_orig[i];//max(tau, x_orig[i]);//x_orig without pivoting, max formulation with pivoting
end for;
......@@ -1530,199 +1510,6 @@ package Common
m[1:nL-1] := Functions.LiquidFunctions.calc_mfromY(Yl);
m[nL] := Yl[nL]/MH2O;
gamma_i := Functions.LiquidFunctions.calc_gamma(T,p,Xl);
gi := Functions.calc_gi(T, p) .* MMX;
gi_id :=P_to_id_*gi;
// gr :=nu_id*gi_id;
gr :=nu_id__*gi;
logK :=-gr/Modelica.Constants.R/T;
a :=gamma_i .* m;
logreacBase :=cat(1,-z[1:ns],log(a)-z[ns+1:ns+nL]);
logreacBase_id :=P_to_id_*logreacBase;
// x_id :=P_to_id_*x;
//isopotential
f[1:nR] :=nu_id__*logreacBase - logK;
// f[1:nR] :=nu_id*logreacBase_id - logK;
//reduced mass balance
// f[nR+1:nF] :=transpose(lambda_mass)*mass - Xred;
// f[nR+1:nF] :=transpose(lambda_nu)*x_id - Xred;//transpose(lambda_mass)*mass - Xred;
f[nR+1:nF] :=transpose(lambda_nu)*x - Xred;//transpose(lambda_mass)*mass - Xred;
end SLE_Tp;
function SLE_Tp_
"SLE initialization according to Leal et al. (2016): calculation of residuals"
input Real[nF] x;
input Real[nX] Xred;
input SI.Temperature T;
input SI.Pressure p;
output Real[nF] f;
protected
SI.MolarEnergy gi[ns+nL];
SI.MolarEnergy gr[nR];
Real[nR] logK;
Real[nL] gamma_i;
Real[nL] m;
Real [nL] a;
SI.MoleFraction[ns+nL] Y;
SI.MassFraction[ns+nL] X;
SI.MoleFraction[nL] Yl;
SI.MassFraction[nL] Xl;
SI.MoleFraction[ns] Ys;
SI.MassFraction[ns] Xs;
Real[ns] slack = 1e-20./x[1:ns];//x[nF+1:nF+ns];
parameter Real eps = 1e-15;
Real[ns+nL] logreacBase;
SI.Mass[ns+nL] mass;
algorithm
//assert(sum(x[1:ns+nL]) > 0, "x is zero in SLE_Tp");
Y :=x[1:ns+nL]/sum(x[1:ns+nL]);
X :=Functions.calc_Xfull(Y);
Ys :=Y[1:ns]/sum(Y[1:ns]);
Yl :=Y[1 + ns:ns + nL]/sum(Y[1 + ns:ns + nL]);
Xs :=Functions.calc_X_M(Ys, MMX[1:ns]);
Xl :=Functions.calc_X_M(Yl, MMX[1 + ns:ns + nL]);
m[1:nL-1] := Functions.LiquidFunctions.calc_mfromY(Yl);
m[nL] := Yl[nL]/MH2O;
gamma_i := Functions.LiquidFunctions.calc_gamma(T,p,Xl);
gi := Functions.calc_gi(T, p) .* MMX;
gr :=nu*gi;
logK :=-gr/Modelica.Constants.R/T;
a :=gamma_i .* m;
logreacBase :=cat(1,-slack,log(a));
mass :=x[1:nF] .* MMX;
//isopotential
f[1:nR] :=nu*logreacBase - logK;
//reduced mass balance
f[nR+1:nF] :=transpose(lambda_mass)*mass - Xred;
end SLE_Tp_;
function SLE_Tp_v1
"SLE initialization according to Leal et al. (2016): calculation of residuals"
input Real[nF] x;
input Real[nX] Xred;
input SI.Temperature T;
input SI.Pressure p;
output Real[nF] f;
protected
SI.MolarEnergy gi[ns+nL];
SI.MolarEnergy gr[nR];
Real[nR] logK;
Real[nL] gamma_i;
Real[nL] m;
Real [nL] a;
SI.MoleFraction[ns+nL] Y;
SI.MassFraction[ns+nL] X;
SI.MoleFraction[nL] Yl;
SI.MassFraction[nL] Xl;
SI.MoleFraction[ns] Ys;
SI.MassFraction[ns] Xs;
Real tau = 1e-30;
SI.MassFraction[ns+nL] z = {if x[i] > 0 then tau/max(tau,x[i]) else 0 for i in 1:ns+nL};//tau./x;
Real[ns+nL] logreacBase;
SI.Mass[ns+nL] mass;
algorithm
assert(sum(x[1:ns+nL]) > 0, "x is zero in SLE_Tp");
Y :=x[1:ns+nL]/sum(x[1:ns+nL]);
X :=Functions.calc_Xfull(Y);
Ys :=Y[1:ns]/sum(Y[1:ns]);
Yl :=Y[1 + ns:ns + nL]/sum(Y[1 + ns:ns + nL]);
Xs :=Functions.calc_X_M(Ys, MMX[1:ns]);
Xl :=Functions.calc_X_M(Yl, MMX[1 + ns:ns + nL]);
m[1:nL-1] := Functions.LiquidFunctions.calc_mfromY(Yl);
m[nL] := Yl[nL]/MH2O;
gamma_i := Functions.LiquidFunctions.calc_gamma(T,p,Xl);
gi := Functions.calc_gi(T, p) .* MMX;
gr :=nu*gi;
logK :=-gr/Modelica.Constants.R/T;
a :=gamma_i .* m;
logreacBase :=cat(1,-z[1:ns],log(a)-z[ns+1:ns+nL]);
mass :=x[1:nF] .* MMX;
//isopotential
f[1:nR] :=nu*logreacBase - logK;
//reduced mass balance
// f[nR+1:nF] :=transpose(lambda_mass)*mass - Xred;
f[nR+1:nF] :=transpose(lambda_)*x - Xred;//transpose(lambda_mass)*mass - Xred;
end SLE_Tp_v1;
function SLE_Tp_v2
"SLE initialization according to Leal et al. (2016): calculation of residuals"
input Real[nF] x;
input Real[nX] Xred;
input SI.Temperature T;
input SI.Pressure p;
output Real[nF] f;
protected
SI.MolarEnergy gi[ns+nL];
SI.MolarEnergy gi_id[ns+nL];
SI.MolarEnergy gr[nR];
Real[nR] logK;
Real[nL] gamma_i;
Real[nL] m;
Real [nL] a;
SI.MoleFraction[ns+nL] Y;
SI.MassFraction[ns+nL] X;
SI.MoleFraction[nL] Yl;
SI.MassFraction[nL] Xl;
SI.MoleFraction[ns] Ys;
SI.MassFraction[ns] Xs;
Real tau = 1e-30;
SI.MassFraction[ns+nL] z = {if x[i] > 0 then tau/max(tau,x[i]) else 0 for i in 1:ns+nL};//tau./x;
Real[ns+nL] logreacBase;
Real[ns+nL] logreacBase_id;
Real[ns+nL] x_id;
SI.Mass[ns+nL] mass;
algorithm
assert(sum(x[1:ns+nL]) > 0, "x is zero in SLE_Tp");
Y :=x[1:ns+nL]/sum(x[1:ns+nL]);
X :=Functions.calc_Xfull(Y);
Ys :=Y[1:ns]/sum(Y[1:ns]);
Yl :=Y[1 + ns:ns + nL]/sum(Y[1 + ns:ns + nL]);
Xs :=Functions.calc_X_M(Ys, MMX[1:ns]);
Xl :=Functions.calc_X_M(Yl, MMX[1 + ns:ns + nL]);
m[1:nL-1] := Functions.LiquidFunctions.calc_mfromY(Yl);
m[nL] := Yl[nL]/MH2O;
gamma_i := Functions.LiquidFunctions.calc_gamma(T,p,Xl);
gi := Functions.calc_gi(T, p) .* MMX;
gi_id :=P_to_id*gi;
......@@ -1735,92 +1522,18 @@ package Common
logreacBase :=cat(1,-z[1:ns],log(a)-z[ns+1:ns+nL]);
logreacBase_id :=P_to_id*logreacBase;
mass :=x[1:nF] .* MMX;
x_id :=P_to_id*x;
// x_id :=P_to_id*x;
//isopotential
f[1:nR] :=nu*logreacBase - logK;
// f[1:nR] :=nu_id*logreacBase_id - logK;
//reduced mass balance
// f[nR+1:nF] :=transpose(lambda_mass)*mass - Xred;
f[nR+1:nF] :=transpose(lambda_id)*x_id - Xred;//transpose(lambda_mass)*mass - Xred;
end SLE_Tp_v2;
function SLE_Tp_v3
"SLE initialization according to Leal et al. (2016): calculation of residuals"
input Real[nF] x;
input Real[nX] Xred;
input SI.Temperature T;
input SI.Pressure p;
output Real[nF] f;
protected
SI.MolarEnergy gi[ns+nL];
SI.MolarEnergy gi_id[ns+nL];
SI.MolarEnergy gr[nR];
Real[nR] logK;
Real[nL] gamma_i;
Real[nL] m;
Real [nL] a;
SI.MoleFraction[ns+nL] Y;
SI.MassFraction[ns+nL] X;
SI.MoleFraction[nL] Yl;
SI.MassFraction[nL] Xl;
SI.MoleFraction[ns] Ys;
SI.MassFraction[ns] Xs;
Real tau = 1e-30;
SI.MassFraction[ns+nL] z;// = {if x[i] > 0 then tau/max(tau,x[i]) else 0 for i in 1:ns+nL};//tau./x;
Real[ns+nL] logreacBase;
Real[ns+nL] logreacBase_id;
Real[ns+nL] x_orig;
SI.Mass[ns+nL] mass;
algorithm
assert(sum(x[1:ns+nL]) > 0, "x is zero in SLE_Tp");
x_orig :=P_to_orig_*x;
z :={if x_orig[i] > 0 then tau/max(tau, x_orig[i]) else 0 for i in 1:ns + nL};
Y :=x_orig[1:ns+nL]/sum(x_orig[1:ns+nL]);
X :=Functions.calc_Xfull(Y);
Ys :=Y[1:ns]/sum(Y[1:ns]);
Yl :=Y[1 + ns:ns + nL]/sum(Y[1 + ns:ns + nL]);
Xs :=Functions.calc_X_M(Ys, MMX[1:ns]);
Xl :=Functions.calc_X_M(Yl, MMX[1 + ns:ns + nL]);
m[1:nL-1] := Functions.LiquidFunctions.calc_mfromY(Yl);
m[nL] := Yl[nL]/MH2O;
gamma_i := Functions.LiquidFunctions.calc_gamma(T,p,Xl);
gi := Functions.calc_gi(T, p) .* MMX;
gi_id :=P_to_id_*gi;
// gr :=nu_id*gi_id;
gr :=nu_id__*gi;
logK :=-gr/Modelica.Constants.R/T;
a :=gamma_i .* m;
logreacBase :=cat(1,-z[1:ns],log(a)-z[ns+1:ns+nL]);
logreacBase_id :=P_to_id_*logreacBase;
// x_id :=P_to_id_*x;
//isopotential
f[1:nR] :=nu_id__*logreacBase - logK;
// f[1:nR] :=nu_id*logreacBase_id - logK;
//reduced mass balance
// f[nR+1:nF] :=transpose(lambda_mass)*mass - Xred;
// f[nR+1:nF] :=transpose(lambda_nu)*x_id - Xred;//transpose(lambda_mass)*mass - Xred;
f[nR+1:nF] :=transpose(lambda_nu)*x - Xred;//transpose(lambda_mass)*mass - Xred;
end SLE_Tp_v3;
f[nR+1:nF] :=transpose(lambda_id)*x - Xred;//transpose(lambda_mass)*mass - Xred;
end SLE_Tp;
end calc_f;
package calc_J
......@@ -1890,7 +1603,7 @@ package Common
Real[ns+nL] x_orig;
algorithm
x_orig :=P_to_orig_*x;
x_orig :=P_to_orig*x;
for i in 1:ns+nL loop
n[i] :=x_orig[i];//max(tau, x_orig[i]);//x_orig without pivoting, max formulation with pivoting
......@@ -1908,7 +1621,7 @@ package Common
diagH[ns+1:ns+nL] :=(1 .- Yl .+ z[ns+1:ns+nL]) ./ n[ns + 1:ns + nL];// + tau ./ n[ns+1:ns+nL] .^ 2;//(1 .- Yl) ./ n[ns + 1:ns + nL] + tau ./ n[ns+1:ns+nL] .^ 2;
H :=diagonal(diagH);//fill(diagH,ns+nL);//
H_id :=transpose(P_to_id_)*H;
H_id :=transpose(P_to_id)*H;
// for i in 1:ns+nL loop
// //isopotential
......@@ -1925,252 +1638,25 @@ package Common
// for r in 1:nR loop
// J[r,:] :=nu[r,:].*diagH;
// end for;
J[1:nR,:] :=nu_id__*H;//nu_id__ .* H;//nu_id__*H;//
J[1:nR,:] :=J[1:nR, :]*transpose(P_to_id_);//transpose(P_nu_id);
J[1:nR,:] :=nu*H;//nu_id__ .* H;//nu_id__*H;//
J[1:nR,:] :=J[1:nR, :]*transpose(P_to_id);//transpose(P_nu_id);
// J[1:nR,:] :=nu_id_ .* H_id;
// J[1:nR,:] :=nu_id .* H_id;
//reduced mass balance
J[nR+1:nF,:] :=transpose(lambda_nu);
J[nR+1:nF,:] :=transpose(lambda_id);
end SLE_Tp;
function SLE_Tp_
"SLE initialization according to Leal et al. (2016): calculation of Jacobian"
input Real[nF] x;
input SI.Temperature T;
input SI.Pressure p;
output Real[nF,nF] J;
protected
SI.MolarEnergy gi[ns+nL];
SI.MoleFraction[ns+nL] Y;
SI.MassFraction[ns+nL] X;
SI.MoleFraction[nL] Yl;
SI.MassFraction[nL] Xl;
SI.MoleFraction[ns] Ys;
SI.MassFraction[ns] Xs;
parameter Real eps = 1e-15;
Real[ns] slack = 1e-20./x[1:ns];//x[nF+1:nF+ns];
Real[ns+nL] logreacBase;
algorithm
Y :=x[1:ns+nL]/sum(x[1:ns+nL]);
X :=Functions.calc_Xfull(Y);
end calc_J;
Ys :=Y[1:ns]/sum(Y[1:ns]);
Yl :=Y[1 + ns:ns + nL]/sum(Y[1 + ns:ns + nL]);
for i in 1:ns+nL loop
//isopotential
for r in 1:nR loop
if i < ns+1 then
J[r,i] := if x[i] > 0 then nu[r,i]*slack[i]/x[i] else 0;
else
J[r,i] :=if nu[r, i] > 0 then nu[r, i]*(1 - Yl[i - ns])/x[i] else 0;
end if;
end for;
//reduced mass balance
for j in 1:ns+nL - nR loop
J[j + nR,i] := lambda_mass[i, j]*MMX[i];
end for;
end for;
end SLE_Tp_;
function SLE_Tp_v1
"SLE initialization according to Leal et al. (2016): calculation of Jacobian"
input Real[nF] x;
input SI.Temperature T;
input SI.Pressure p;
output Real[nF,nF] J;
protected
SI.AmountOfSubstance[ns+nL] n;
SI.MoleFraction[ns+nL] Y;
SI.MassFraction[ns+nL] X;
SI.MoleFraction[nL] Yl;
SI.MassFraction[nL] Xl;
SI.MoleFraction[ns] Ys;
SI.MassFraction[ns] Xs;
Real tau = 1e-20;
SI.MassFraction[ns+nL] z;
Real[ns+nL,ns+nL] H;
Real[ns+nL] diagH;
algorithm
Y :=x[1:ns+nL]/sum(x[1:ns+nL]);
X :=Functions.calc_Xfull(Y);
Ys :=Y[1:ns]/sum(Y[1:ns]);
Yl :=Y[1 + ns:ns + nL]/sum(Y[1 + ns:ns + nL]);
for i in 1:ns+nL loop
n[i] :=max(tau, x[i]);
end for;
z :=tau./n;
diagH[1:ns] :=tau ./ n[1:ns] .^ 2;
diagH[ns+1:ns+nL] :=(1 .- Yl) ./ n[ns + 1:ns + nL] + tau ./ n[ns+1:ns+nL] .^ 2;
H :=fill(diagH,ns+nL);
// for i in 1:ns+nL loop
// //isopotential
// for r in 1:nR loop
// if i < ns+1 then
// J[r,i] := nu[r,i]*tau/x[i]^2 else 0;//nu[r,i]*z[i]/n[i] else 0;//if abs(nu[r,i]) > 0 then nu[r,i]*z[i]/n[i] else 0;
// else
// J[r,i] := if abs(nu[r,i])>0 then nu[r, i]*(1 - Yl[i - ns])/x[i] + nu[r,i]*tau/x[i]^2 else 0;//nu[r, i]*(1 - Yl[i - ns])/n[i] + nu[r,i]*z[i]/n[i] else 0;//if abs(nu[r,i]) > 0 then nu[r, i]*(1 - Yl[i - ns])/n[i] + nu[r,i]*z[i]/n[i] else 0;
// end if;
// end for;
// end for;
//isopotential
// for r in 1:nR loop
// J[r,:] :=nu[r,:].*diagH;
// end for;
J[1:nR,:] :=nu .* H;
//reduced mass balance
J[nR+1:nF,:] :=transpose(lambda_);
end SLE_Tp_v1;
function SLE_Tp_v2
"SLE initialization according to Leal et al. (2016): calculation of Jacobian"