Update on LiquidPhase initialization according to Leal et al. (2016)

ElectrolyteMedia/Media/LiquidPhase/Common/MixtureLiquid/Initialization/calc_J/LE_Tp.mo

@@ -7,34 +7,39 @@ function LE_Tp
   output Real[nL,nL] J;
 
 protected
+  SI.AmountOfSubstance[nF] n;
   SI.MoleFraction[nL] Y;
   SI.MassFraction[nL] X;
   Real[nL] gamma_i;
   Real[nL] m;
+  Real tau = 1e-37;
+  SI.MassFraction[nF] z;
+  Real[nF,nF] H;
+  Real[nF,nF] H_id;
+  Real[nF] diagH;
+  Real[nF] x_orig;
 algorithm
-  assert(sum(x)>0,"Input error: x is zero");
-  Y :=x/sum(x);
 
-  m[1:nL-1] := Functions.calc_mfromY(Y);
-  m[nL] := Y[nL]/MH2O;
+  x_orig :=P_to_orig*x;
 
+  for i in 1:nF loop
+    n[i] :=x_orig[i];//max(tau, x_orig[i]);//x_orig without pivoting, max formulation with pivoting
+  end for;
+
+  Y :=n[1:nF]/sum(n[1:nF]);
   X :=Functions.calc_X(Y);
 
-  gamma_i := Functions.calc_gamma(T,p,X);
-
-  for i in  1:nL loop
-    //isopotential
-    for r in 1:nR loop
-      if i < nL then
-        J[r, i] := if abs(nu[r,i]) > 0 then nu[r,i]/x[i] - nu[r,nL]/sum(x) else - nu[r,nL]/sum(x);
-      else
-        J[r,i] :=1/x[i]*(nu[r,i]-sum(nu[r,1:nL-1])) - nu[r,i]/sum(x);
-      end if;
-    end for;
-
-    //reduced mass balance
-    for j in 1:nL - nR loop
-      J[j + nR, i] := lambda_mass[i, j]*MMX[i];
-    end for;
-  end for;
+  z :=tau./n;
+
+  diagH :=(1 .- Y .+ z) ./ n;
+
+  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 LE_Tp;

ElectrolyteMedia/Media/LiquidPhase/Common/MixtureLiquid/Initialization/calc_f/LE_Tp.mo

@@ -8,19 +8,26 @@ function LE_Tp
   output Real[nL] f;
 
 protected
-  SI.MolarEnergy gi[nL];
+  SI.MolarEnergy gi[nF];
   SI.MolarEnergy gr[nR];
   Real[nR] logK;
   Real[nL] gamma_i;
   Real prod;
   Real[nL] m;
   Real[nL] loga;
+  Real tau = 1e-37;
+  Real[nF] z;
   SI.MoleFraction[nL] Y;
   SI.MassFraction[nL] X;
-  SI.Mass[nL] mass;
+  Real[nF] x_orig "moles in original order";
 algorithm
 
-  Y :=x/sum(x);
+  assert(sum(x[1:nF]) > 0, "x is zero in LE_Tp");
+
+  x_orig :=P_to_orig*x;
+  z :=tau ./ x_orig;
+
+  Y :=x_orig/sum(x_orig);
 
   m[1:nL-1] := Functions.calc_mfromY(Y);
   m[nL] := Y[nL]/MH2O;
@@ -35,11 +42,10 @@ algorithm
   logK :=-gr/Modelica.Constants.R/T;
 
   loga :=log(gamma_i .* m);
-  mass :=x .* MMX;
 
   //isopotential
-  f[1:nR] :=nu*loga - logK;
+  f[1:nR] :=nu*(loga - z) - logK;
 
   //reduced mass balance
-  f[nR+1:nF] :=transpose(lambda_mass)*mass - Xred;
+  f[nR+1:nF] :=transpose(lambda_id)*x - Xred;
 end LE_Tp;

ElectrolyteMedia/Media/LiquidPhase/Common/MixtureLiquid/X_pTXred.mo

@@ -4,90 +4,79 @@ function X_pTXred
   input SI.Pressure p;
   input SI.Temperature T;
   input MassFraction[nX] Xred;
-  input MassFraction[nL] Xinit = zeros(nL);
-  output Real[nL] Xfull;
+  input MassFraction[nF] Xinit = zeros(nF);
+  output Real[nF] x_;
 
 protected
   parameter Real init = 1e-10;
   parameter Real eps = 1e-6;
-  parameter Integer[nR] firstnonzero = {Modelica.Math.BooleanVectors.firstTrueIndex({nu[r,i] <> 0 for i in 1:nL}) for r in 1:nR};
-  parameter Integer[nR] secondnonzero = {Modelica.Math.BooleanVectors.firstTrueIndex({if i>firstnonzero[r] then nu[r,i] <> 0 else false for i in 1:nL}) for r in 1:nR};
 
   MassFraction[nX] Xrednonzero;
+  MassFraction[nX] Xred_;
 
-  Real[nL,nL] J;
-  Real[nL] f;
-  Real[nL] x1;
-  Real[nL] x2;
+  Real[nF,nF] J;
+  Real[nF] f;
+  Real[nF] x1;
+  Real[nF] x2;
 
-  Real[nR,nL] B;
-  Real[nR] b;
+  SI.MoleFraction[nF] Yfull;
+  SI.MassFraction[nF] Xfull;
 
-  SI.MassFraction[nL] Xfull_in;
-  SI.Mass[nL] mfull_in;
-  SI.MoleFraction[nL] Yfull_in;
-  SI.AmountOfSubstance[nL] nfull_in;
-
-  SI.MoleFraction[nL] Yfull;
-
-  SI.MoleFraction[nL] Yinit;
-  SI.AmountOfSubstance[nL] ninit;
-  SI.Mass[nL] minit;
+  SI.MoleFraction[nF] Yinit;
+  SI.AmountOfSubstance[nF] ninit;
+  SI.Mass[nF] minit;
   Real scale;
 
   Boolean solutionfound = false;
 
+  Real[nX] Xs;
 algorithm
-
   //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;
 
+  //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;
+
+  Xred_ :=transpose(lambda_mass_id_orig)*Xfull;
+
   // find initial guess x1 for Newton solver
-  for r in 1:nR loop
-    B[r,secondnonzero[r]] :=1;
-    b[r] :=0;
-  end for;
-  Xfull_in := ElectrolyteMedia.Media.Common.equalityLeastSquares(
-    transpose(lambda_mass),
-    Xred,
-    B,
-    b);
-  // nfull_in with positive elements only to fulfill transpose(lambda_mass)*mfull_in = Xred;
-  if sum(Xfull_in) > 0 then
-    mfull_in :=Xfull_in;
-    nfull_in :=mfull_in ./ MMX;
-    for i in 1:nF loop
-      x1[i] :=max(init, nfull_in[i]);
-    end for;
-  elseif sum(Xinit)> 0 then
+  if sum(Xinit)> 0 then
   // ninit from initial guess provided as input (Xinit)
-    scale :=sum(transpose(lambda_mass)*Xinit);
-    minit :=Xinit*scale;
-    Yinit :=Functions.calc_Y(Xinit);
-    ninit :=minit./MMX;
-    for i in 1:nF loop
-      x1[i] :=max(init, ninit[i]);
-    end for;
+     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
-    x1 :=fill(1/nF, nF);
+    for i in 1:nF loop
+      x1[i] :=max(init, Xfull[i]/MMX[i]);
+    end for;
   end if;
 
   // Newton algorithm
   solutionfound :=false;
+  x1 :=P_to_id*x1;
   while not solutionfound loop
-    f :=Initialization.calc_f.LE_Tp(x1,Xrednonzero,T,p);
+    f :=Initialization.calc_f.LE_Tp(x1,Xred_,T,p);
     J :=Initialization.calc_J.LE_Tp(x1,T,p);
-    x2 := Initialization.NewtonStep(x1,f,J,nL);
+    x2 := Initialization.NewtonStep(x1,f,J,nF);
+
+    // check convergence
     if Modelica.Math.Vectors.norm(f,Modelica.Constants.inf) < eps then
-      solutionfound :=true;
+      solutionfound:=true;
     end if;
     x1 :=x2;
   end while;
+  x2 :=P_to_orig*x2;
+  Yfull :=x2[1:nF]/sum(x2[1:nF]);
+  Xfull :=moleToMassFractions(Yfull, MMX);
 
-  Yfull :=x2/sum(x2);
-  Xfull :=Functions.calc_X(Yfull);
+  x_ :=Xfull;
 
 end X_pTXred;

ElectrolyteMedia/Media/LiquidPhase/Common/MixtureLiquid/package.mo

@@ -36,8 +36,18 @@ package MixtureLiquid "Medium model of a mixture of liquids based on Modelica li
   constant Real[nR,nF] nu = userInterface.nu "Stoichiometry matrix of gas-liquid and dissociation equilibria from user interface";
   constant Real[nR,nF] nu_mass = Functions.Reaction.calc_nu_mass(nu) "Mass based stoichiometry matrix of gas-liquid and dissociation equilibria from user interface";
 
-  constant Real[:,:] lambda= Media.Common.Reaction.calc_lambda(nu, nR, nF) "Null space of stoichiometry matrix";
-  constant Real[:,:] lambda_mass= {{lambda[i,j]/MMX[i] for j in 1:nX} for i in 1:nL}*1000 "Null space of mass based stoichiometry matrix";
+//   constant Real[:,:] lambda= Media.Common.Reaction.calc_lambda(nu, nR, nF) "Null space of stoichiometry matrix";
+//   constant Real[:,:] lambda_mass= {{lambda[i,j]/MMX[i] for j in 1:nX} for i in 1:nL}*1000 "Null space of mass based stoichiometry matrix";
+  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}/1000 "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 "Gas constant of H2O";

ElectrolyteMedia/Media/LiquidPhase/Common/MixtureLiquid/package.order

@@ -13,6 +13,13 @@ nu
 nu_mass
 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

ElectrolyteMedia/Media/SolidLiquidPhase/Common/MixtureLiquid/X_pTXred.mo

@@ -28,8 +28,6 @@ protected
   SI.MoleFraction[nF] Yfull;
   SI.MassFraction[nF] Xfull;
 
-  Integer index;
-
   Real[nX] Xs;
 algorithm
   //ensure no species amount to be zero for numerical issues

ElectrolyteMedia/Tests/LiquidPhase/NaOH_HCL_titration.mo

@@ -31,5 +31,5 @@ model NaOH_HCL_titration "Implementation of partial_Titration"
     mDot_titrant=2.5e-5) annotation (Placement(transformation(extent={{-10,-14},{10,6}})));
    annotation (Icon(coordinateSystem(preserveAspectRatio=false)), Diagram(
         coordinateSystem(preserveAspectRatio=false)),
-    experiment(StopTime=2000, __Dymola_NumberOfIntervals=4000));
+    experiment(StopTime=2000, Tolerance=1e-06));
 end NaOH_HCL_titration;