Commit 134ee770 authored by Andreas Michael Bremen's avatar Andreas Michael Bremen
Browse files

LiquidPhase initialization according to Leal et al. (2016)

Minor changes to SolidLiquidPhase
parent 9a393c17
within ElectrolyteMedia.Media.GasLiquidPhase.Common.Functions;
function calc_Xfull
"calculates full mass fraction vector from full mole fraction vector"
input SI.MoleFraction[nLfunfun+nGfunfun] Y;
output SI.MassFraction [nLfunfun+nGfunfun] X;
protected
Real[nLfunfun+nGfunfun] Y_(unit="mol/kg");
algorithm
for i in nGfunfun+1:nGfunfun+nLfunfun-1 loop
Y_[i] :=Y[i]*dataLfunfun[i-nGfunfun].MM;
end for;
for i in 1:nGfunfun loop
Y_[i] :=Y[i]* dataGfunfun[i].MM;//0.05844;//
end for;
Y_[nGfunfun+nLfunfun] :=Y[nGfunfun+nLfunfun]*IF97.MH2O;
X :=Y_/sum(Y_);
end calc_Xfull;
......@@ -30,3 +30,4 @@ LiquidFunctions
GasFunctions
moleToMassFractions
massToMoleFractions
calc_Xfull
within ElectrolyteMedia.Media.GasLiquidPhase.Common.MixtureLiquid.Initialization.calc_J;
function GLE_Tp
"Gas-liquid and dissociation equilibrium with constraints on inert solute molalities and moles of gas phase"
input Real[nG+nL] x "mole fractions of gas and liquid phase for each phase separately";
input Real[nG+nL] inert "mole fraction of inert gas phase species and molalities of inert liquid phase species";
"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[nG + nL,nG + nL] J;
output Real[nF,nF] J;
protected
Boolean[nG+nL] isinert = {sum(abs(nu[:,i])) < Modelica.Constants.eps for i in 1:nG+nL};
Integer index;
SI.AmountOfSubstance[nG+nL] n;
SI.MoleFraction[nG+nL] Y;
SI.MassFraction[nG+nL] X;
SI.MoleFraction[nL] Yl;
SI.MassFraction[nL] Xl;
SI.MoleFraction[nG] Yg;
SI.MassFraction[nG] Xg;
Real tau = 1e-37;
SI.MassFraction[nG+nL] z;
Real[nG+nL,nG+nL] H;
Real[nG+nL,nG+nL] H_id;
Real[nG+nL] diagH;
Real[nG+nL] x_orig;
algorithm
for i in 1:nG + nL loop
//isopotential
for r in 1:nR loop
if i < nG+1 then
J[r,i] := if nu[r,i]<> 0 then nu[r, i]/x[i] else 0;
elseif i == nG+nL then
J[r,i] :=-sum({nu[r, i_] for i_ in 1 + nG:nG + nL - 1})/x[i] + nu[r, i]/x[i];
else J[r,i] :=if abs(nu[r,i]) > 0 then nu[r, i]/x[i] else 0;
end if;
end for;
//nR+1: closure condition gas phase
//nR+2: closure condition liquid phase
//nR+3: charge balance liquid phase
if i < nG+1 then
J[nR+1,i] :=1;
J[nR+2,i] :=0;
J[nR+3,i] :=0;
elseif i == nG+nL then
J[nR+1,i] :=0;
J[nR+2,i] :=1;
J[nR+3,i] :=0;
else
J[nR+1,i] :=0;
J[nR+2,i] :=1;
J[nR+3,i] :=datal[i-nG].z;
end if;
end for;
//constraints on inert solute molalities and mole fraction of gas phase
index :=4;
for i in 1:nG+nL-1 loop
if isinert[i] then
if i < nG+1 then
J[nR+index,i] :=1000;
else
J[nR+index,i] :=1;
J[nR+index,nG+nL] :=-MH2O*inert[i];
end if;
index :=index + 1;
if index > nG+nL-nR then
break;
end if;
end if;
x_orig :=P_to_orig*x;
for i in 1:nG+nL loop
n[i] :=x_orig[i];//max(tau, x_orig[i]);//x_orig without pivoting, max formulation with pivoting
end for;
Y :=n[1:nG+nL]/sum(n[1:nG+nL]);
X :=Functions.calc_Xfull(Y);
Yg :=Y[1:nG]/sum(Y[1:nG]);
Yl :=Y[1 + nG:nG + nL]/sum(Y[1 + nG:nG + nL]);
z :=tau./n;
diagH[1:nG] :=(1 .- Yg .+ z[1:nG]) ./ n[1:nG];
diagH[nG+1:nG+nL] :=(1 .- Yl .+ z[nG+1:nG+nL]) ./ n[nG + 1:nG + nL];
H :=diagonal(diagH);
H_id :=transpose(P_to_id)*H;
//isopotential
J[1:nR,:] :=nu*H;
J[1:nR,:] :=J[1:nR, :]*transpose(P_to_id);
//reduced mass balance
J[nR+1:nF,:] :=transpose(lambda_id);
end GLE_Tp;
within ElectrolyteMedia.Media.GasLiquidPhase.Common.MixtureLiquid.Initialization.calc_f;
function GLE_Tp
"Gas-liquid and dissociation equilibrium with constraints on solute molalities and moles of gas phase"
input Real[nG+nL] x "mole fractions of gas and liquid phase for each phase separately";
input Real[nG+nL] inert "mole fraction of inert gas phase species and molalities of inert liquid phase species";
"GLE initialization according to Leal et al. (2016): calculation of residuals"
input Real[nG+nL] x;
input Real[nX] Xred;
input SI.Temperature T;
input SI.Pressure p;
output Real[nG+nL] f;
protected
SI.SpecificEnergy gi[nG+nL];
SI.MolarEnergy gi[nG+nL];
SI.MolarEnergy gim[nG+nL];
SI.MolarEnergy gr[nR];
Real[nR] K;
Real[nR] logK;
Real[nL] gamma_i;
Real[nL] m;
SI.MoleFraction[nG] Yg = x[1:nG]/sum(x[1:nG]);
SI.MoleFraction[nL] Yl = x[1+nG:nG+nL]/sum(x[1+nG:nG+nL]);
SI.MassFraction[nG] Xg = Functions.calc_X_M(Yg, MMX[1:nG]);
SI.MassFraction[nL] Xl = Functions.LiquidFunctions.calc_X(Yl);
Real [nL] a;
SI.MoleFraction[nG] Yg;
SI.MoleFraction[nL] Yl;
SI.MassFraction[nG] Xg;
SI.MassFraction[nL] Xl;
Real[nG] phi;
Integer index = 1;
Real[nG] fug;
SI.Density dg;
Real tau = 1e-37;
SI.MassFraction[nG+nL] z;
Real[nG+nL] logreacBase;
Real[nG+nL] x_orig "moles in original order";
algorithm
assert(sum(x) > 0, "x is zero in GLE_Tp");
gi :=Functions.calc_gi(T, p);
for i in 1:nG+nL loop
gim[i] :=gi[i]*MMX[i];
end for;
gamma_i := Functions.LiquidFunctions.calc_gamma(
T,
p,
Xl);
dg :=dg_TpX(
T,
p,
Xg);
x_orig :=P_to_orig*x;
z :=tau ./ x_orig;
Yg :=x_orig[1:nG]/sum(x_orig[1:nG]);
Yl :=x_orig[1 + nG:nG + nL]/sum(x_orig[1 + nG:nG + nL]);
Xg :=Functions.calc_X_M(Yg, MMX[1:nG]);
Xl :=Functions.LiquidFunctions.calc_X(Yl);
gamma_i := Functions.LiquidFunctions.calc_gamma(T,p,Xl);
dg :=dg_TpX(T,p,Xg);
phi :=Functions.GasFunctions.calc_phi(T,dg,Xg);
m[1:nL-1] :=Functions.LiquidFunctions.calc_mfromY(x[1+ nG:nG + nL]);
m[nL] :=x[nG+nL]/MH2O;
m[1:nL-1] :=Functions.LiquidFunctions.calc_mfromY(x_orig[1+ nG:nG + nL]);
m[nL] :=Xl[nL]/MH2O;
for j in 1:nR loop
gr[j] :=sum({nu[j, i]*gim[i] for i in 1:nG+nL});
K[j] :=exp(-gr[j]/Modelica.Constants.R/T);
end for;
gi :=Functions.calc_gi(T, p) .* MMX;
gr :=nu*gi;
logK :=-gr/Modelica.Constants.R/T;
//isopotential
for r in 1:nR loop
f[r] :=(sum({if nu[r,i_] <> 0 then log(x[i_]*phi[i_]*p/prefg)*nu[r,i_] else 0 for i_ in 1:nG})+sum({if nu[r,i_+nG] <> 0 then log(gamma_i[i_]*m[i_])*nu[r, i_+nG] else 0 for i_ in 1:nL}) - log(K[r]));
end for;
a :=gamma_i .* m;
fug :=Yg .* phi*p/prefg;
//nR+1: closure condition gas phase
f[nR+1] :=sum(x[1:nG]) - 1;
logreacBase :=cat(1,log(fug) - z[1:nG], log(a) - z[nG + 1:nG + nL]);
//nR+2: closure condition liquid phase
f[nR+2] :=sum(x[1 + nG:nG + nL]) - 1;
//isopotential
f[1:nR] :=nu*logreacBase - logK;
//nR+3: charge balance liquid phase
f[nR+3] :=datal[:].z*x[1 + nG:nG + nL - 1];
//reduced mass balance
f[nR+1:nF] :=transpose(lambda_id)*x - Xred;
//constraints on inert solute molalities and mole fraction of gas phase
index :=4;
for i in 1:nG+nL loop
if sum(abs(nu[:,i])) < Modelica.Constants.eps then
if i < nG+1 then
f[nR+index] :=1000*(x[i] - inert[i]);
else
f[nR+index] :=x[i] - x[nG + nL]*MH2O*inert[i];
end if;
index :=index + 1;
if index > nG+nL-nR then
break;
end if;
end if;
end for;
// //isopotential
// for r in 1:nR loop
// f[r] :=(sum({if nu[r,i_] <> 0 then log(x[i_]*phi[i_]*p/prefg)*nu[r,i_] else 0 for i_ in 1:nG})+sum({if nu[r,i_+nG] <> 0 then log(gamma_i[i_]*m[i_])*nu[r, i_+nG] else 0 for i_ in 1:nL}) - log(K[r]));
// end for;
//
// //nR+1: closure condition gas phase
// f[nR+1] :=sum(x[1:nG]) - 1;
//
// //nR+2: closure condition liquid phase
// f[nR+2] :=sum(x[1 + nG:nG + nL]) - 1;
//
// //nR+3: charge balance liquid phase
// f[nR+3] :=datal[:].z*x[1 + nG:nG + nL - 1];
//
// //constraints on inert solute molalities and mole fraction of gas phase
// index :=4;
// for i in 1:nG+nL loop
// if sum(abs(nu[:,i])) < Modelica.Constants.eps then
// if i < nG+1 then
// f[nR+index] :=1000*(x[i] - inert[i]);
// else
// f[nR+index] :=x[i] - x[nG + nL]*MH2O*inert[i];
// end if;
// index :=index + 1;
// if index > nG+nL-nR then
// break;
// end if;
// end if;
// end for;
end GLE_Tp;
......@@ -32,11 +32,11 @@ record UserInterface
//Reaction
parameter Integer nR_GL = 1 "Number of gas-liquid equilibria" annotation(Dialog(tab="Reaction"));
Real[nR_GL,nG+nL] nu_GL = zeros(nR_GL,nG+nL) "Stoichiometry matrix of gas-liquid and dissociation equilibria" annotation(Dialog(tab="Reaction"));
Real[nR_GL,nG+nL] nu_GL = cat(2,identity(nR_GL),-ones(nR_GL,nG+nL-nR_GL)) "Stoichiometry matrix of gas-liquid and dissociation equilibria" annotation(Dialog(tab="Reaction"));
parameter Integer nR_L = 1 "Number of dissociation equilibria" annotation(Dialog(tab="Reaction"));
Real[nR_L,nL] nu_L = zeros(nR_L,nL) "Stoichiometry matrix of dissociation equilibria" annotation(Dialog(tab="Reaction"));
Real[nR_L,nL] nu_L = cat(2,identity(nR_L),-ones(nR_L,nL-nR_L)) "Stoichiometry matrix of dissociation equilibria" annotation(Dialog(tab="Reaction"));
final parameter Integer nR = nR_GL + nR_L "Number of gas-liquid and dissociation equilibria";
......
......@@ -8,7 +8,11 @@ function X_pTXred
output MassFraction[nF] x_;
protected
parameter Real eps = 1e-5;
parameter Real init = 1e-20;
parameter Real eps = 1e-8;
MassFraction[nX] Xrednonzero;
MassFraction[nX] Xred_;
Real[nF,nF] J;
Real[nF] f;
......@@ -17,142 +21,61 @@ protected
Boolean solutionfound = false;
parameter Integer[nR] firstnonzero = {Modelica.Math.BooleanVectors.firstTrueIndex({nu[r,i] <> 0 for i in 1:nF}) for r in 1:nR};
Boolean[nF] isinert = {sum(abs(nu[:,i])) < Modelica.Constants.eps for i in 1:nG+nL};
Real lfrac;
Real[nF] inert;
Real[nR,nF] B;
Real[nR] b;
SI.MassFraction[nG] Xg;
SI.MoleFraction[nG] Yg;
SI.MassFraction[nL] Xl;
SI.MoleFraction[nL] Yl;
Real[nL-1] ml;
SI.MassFraction[nF] Xfull_in;
SI.AmountOfSubstance[nF] Xnotinert;
SI.AmountOfSubstance[nX] Xrednotinert;
Real[nX] Xred_;
SI.AmountOfSubstance[nF] ninit;
SI.Mass[nF] minit;
Real scale;
SI.MoleFraction[nF] Yfull;
SI.MassFraction[nF] Xfull;
SI.MoleFraction[nF] Yinit;
Integer index;
Real[nX] Xs;
algorithm
// find initial guess nfull_in with positive elements only and 1kg in total
for r in 1:nR loop
B[r,firstnonzero[r]] :=1;
b[r] :=0;
//ensure no species amount to be zero for numerical issues
for i in 1:nX loop
Xrednonzero[i] :=max(Modelica.Constants.eps, Xred[i]);
end for;
Xfull_in :=Modelica.Math.Matrices.equalityLeastSquares(transpose(lambda_mass), Xred,B,b);
Xfull_in :=Xfull_in/sum(Xfull_in);
// initial guess for Newton solver
if sum(Xinit)> 0 then
Yinit :=Functions.calc_Y(Xinit);
Yinit[1:nG] :=Yinit[1:nG]/sum(Yinit[1:nG]);
Yinit[1+nG:nF] :=Yinit[1 + nG:nF]/sum(Yinit[1 + nG:nF]);
end if;
//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;
for i in 1:nG+nL loop
if not isinert[i] then
Xnotinert[i] :=Xfull_in[i];
end if;
end for;
Xred_ :=transpose(lambda_mass_id_orig)*Xfull;
//find index for calculation of liquid fraction with non-inert species
Xrednotinert :=transpose(lambda_mass)*Xnotinert;
if index < Modelica.Constants.eps then
index :=Modelica.Math.BooleanVectors.firstTrueIndex({Xrednotinert[i] > 0 for i in 1:nX});
// 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_id_orig)*Xinit);
minit :=Xinit*scale;
ninit :=minit./MMX;
for i in 1:nF loop
x1[i] :=max(init, ninit[i]);
end for;
else
// else generic initial guess
for i in 1:nF loop
x1[i] :=max(init, Xfull[i]/MMX[i]);
end for;
end if;
Xred_ :=transpose(lambda_mass)*Xfull_in;
// specify inert: fixed mole fractions of inert gas species and fixed molalities of inert liquid species
Xg :=Xfull_in[1:nG]/(sum(Xfull_in[1:nG]) + Modelica.Constants.eps);
Yg :=if sum(Xg) > 0 then Functions.GasFunctions.calc_Y(Xg) else zeros(nG);
Xl :=Xfull_in[1+nG:nF]/sum(Xfull_in[1+nG:nF]);
ml :=Functions.LiquidFunctions.calc_mfromX(Xl);
for i in 1:nG+nL loop
if isinert[i] then
if i < nG+1 then
inert[i] :=min(0.9,Yg[i]);// leads to errors if Yg[i] = 1
else
inert[i] :=ml[i-nG];
end if;
end if;
end for;
// Newton algorithm
solutionfound :=false;
x1 :=P_to_id*x1;
while not solutionfound loop
// calculate mole fractions of gas phase and liquid phase
Yfull :=X_pTinert(p,T,inert,Yinit);
f :=Initialization.calc_f.GLE_Tp(x1,Xred_,T,p);
J :=Initialization.calc_J.GLE_Tp(x1,T,p);
x2 := Initialization.NewtonStep(x1,f,J,nF);
// determine liquid fraction lfrac from noninert species
x1[1:nG] :=Functions.GasFunctions.calc_X(Yfull[1:nG]);
x1[1+nG:nG+nL] :=Functions.LiquidFunctions.calc_X(Yfull[1+nG:nF]);
lfrac :=(Xred_[index] - x1[1:nG]*lambda_mass[1:nG, index])/(x1[1 + nG:nG + nL]*
lambda_mass[1 + nG:nG + nL, index] - x1[1:nG]*lambda_mass[1:nG, index]);
for i in 1:nG+nL loop
if i < nG+1 then
x1[i] :=x1[i]*(1 - lfrac);
else
x1[i] :=x1[ i]*lfrac;
end if;
end for;
// lfrac needs to be between 0 and 1
assert(lfrac < 1, "Input mass fraction is not in two phase region --> only liquid phase!",AssertionLevel.error);
assert(lfrac > 0, "Input mass fraction is not in two phase region --> only gas phase!",AssertionLevel.error);
// define break criterion
if Modelica.Math.Vectors.norm((x2-x1)./(x2),Modelica.Constants.inf)< eps then
solutionfound :=true;
// check convergence
if Modelica.Math.Vectors.norm(f,Modelica.Constants.inf) < eps then
solutionfound:=true;
end if;
//store iteration result in x2
x2 :=x1;
// update inert species for next iteration (mole fractions for gaseous species and molalities for liquid species)
Xg :=x1[1:nG];
for i in 1:nG loop
if isinert[i] then
Xg[i] :=Xfull_in[i];
end if;
end for;
Xg :=Xg/sum(Xg);
Yg :=Functions.GasFunctions.calc_Y(Xg);
Xl :=x1[1 + nG:nF];
Yl :=Functions.LiquidFunctions.calc_Y(Xl);
for i in 1:nL loop
if isinert[i+nG] then
Xl[i] :=Xfull_in[i + nG];
end if;
end for;
Xl :=Xl/sum(Xl);
ml :=Functions.LiquidFunctions.calc_mfromX(Xl);
// fixed mole fractions of inert gas species and fixed molalities of inert liquid species
for i in 1:nG+nL-1 loop
if isinert[i] then
if i < nG+1 then
inert[i] :=Yg[i];
else
inert[i] :=ml[i-nG];
end if;
end if;
end for;
//provide start value for next iteration
Yinit :=cat(1,Yg,Yl);
x1 :=x2;
end while;
x2 :=P_to_orig*x2;
Yfull :=x2[1:nF]/sum(x2[1:nF]);
Xfull :=moleToMassFractions(Yfull, MMX);
x_ :=x1;
x_ :=Xfull;
end X_pTXred;
within ElectrolyteMedia.Media.GasLiquidPhase.Common.MixtureLiquid;
function X_pTinert
"Gas-liquid and liquid dissociation equilibrium with output mole fractions of both liquid and gas phase separately. Fixed mole fractions of inert gas species and fixed molalities of inert liquid species"
input SI.Pressure p;
input SI.Temperature T;
input Real[nG+nL] inert;
input MoleFraction[nG+nL] xinit = zeros(nF);
output MoleFraction[nG+nL] x_ "Mass fraction for each gas and liquid phase";
protected
Real[nG + nL,nG + nL] J;
Real[nG + nL] f;
Real[2,nG + nL] x = ones(2,nG + nL);
Real eps = 1e-4;
Boolean solutionfound = false;
algorithm
if sum(xinit) > 0 then
x[1,:] :=xinit;
else
for i in 1:nG loop
x[1,i] :=1/nG;
end for;
x[1,1+nG:nG+nL-1] :=fill(1e-7,nL-1);
x[1,nG+nL] :=1;
end if;
while not solutionfound loop
f :=Initialization.calc_f.GLE_Tp(
x[1, :],
inert,
T,
p);
J :=Initialization.calc_J.GLE_Tp(
x[1, :],
inert,
T,
p);
x[2, :] :=Initialization.NewtonStep(
x[1, :],
f,
J,
nG + nL);
if Modelica.Math.Vectors.norm((x[2,:]-x[1,:])./(x[2,:]),Modelica.Constants.inf)< eps then
solutionfound :=true;
end if;
x[1,:] :=x[2,:];
end while;
x_ :=x[1, :];
end X_pTinert;
......@@ -46,12 +46,16 @@ constant Media.GasLiquidPhase.Common.GasInteractionDataRecord
constant UserInterface userInterface annotation (Placement(transformation(extent={{-10,-8},{10,12}})));
constant Real[:,:] lambda=
Media.Common.Reaction.calc_lambda(
nu,
nR,
nF);
constant Real[:,:] lambda_mass={{lambda[i,j]/MMX[i] for j in 1:nX} for i in 1:nF}*1000 "Null space of mass based stoichiometry matrix";
constant Real[:,:] lambda= Media.Common.Reaction.calc_lambda(nu,nR,nF);
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";
constant Real[nF,nX] lambda_mass_id=P_to_id*lambda_mass "Transformed nullspace of mass based stoichiometry matrix";
constant Real[nF,nX] lambda_mass_id_orig = P_to_orig * {{lambda_id[i,j]/MMX_id[i] for j in 1:nX} for i in 1:nF} "Mass based nullspace back transformed from identity transofrmation";
constant Real[nF] MMX_id = P_to_id*MMX;
constant MolarMass MH2O = IF97.MH2O "Molar mass of solvent";
constant SpecificHeatCapacity RH2O = IF97.RH2O;
......
......@@ -17,6 +17,13 @@ nu_mass
userInterface
lambda
lambda_mass
nu_id
P_to_orig
P_to_id
lambda_id
lambda_mass_id
lambda_mass_id_orig
MMX_id
MH2O
RH2O
MMX
......@@ -47,7 +54,6 @@ T_phX
T_psX
p_TdX
X_pTXred
X_pTinert
XT_phXred
XT_psXred
Xp_dTXred
......
......@@ -12,7 +12,7 @@ protected
Real[nNS] alpha;
Real delta = 0.9999;
Real[nNS,nNS] invJ= Modelica.Math.Matrices.inv(J);
Real[nNS] delta_x = invJ*f;//Modelica.Math.Matrices.solve(J,f);//
Real[nNS] delta_x = invJ*f;
algorithm
for i in 1:nNS loop
......@@ -24,12 +24,4 @@ algorithm
end for;
beta :=min(alpha);
x_ :=x - beta*delta_x;
// while min(x_temp) < 0 loop
// x_temp :=x - 1/beta*invJ*f;
// beta :=beta *10;
// end while;
// assert(min(x_temp)>0,"Newton step too small");
// x_ :=x_temp;
</