Commit c5623947 authored by Andreas Michael Bremen's avatar Andreas Michael Bremen
Browse files

further cleanup of SolidLiquidPhase and resolved bug from previous commit

parent 973c892b
within ElectrolyteMedia.Media.Common.Reaction;
function calc_nu_mass "Mass based stoichiometry"
input Real[:,:] nu;
input Real[size(nu,2)] MMX;
output Real[size(nu,1),size(nu,2)] nu_mass;
algorithm
for i in 1:size(nu,2) loop
nu_mass[:, i] := nu[:, i]*MMX[i];
end for;
end calc_nu_mass;
......@@ -5,3 +5,4 @@ calc_L
calc_U
calc_P_nu_id
calc_lambda_id
calc_nu_mass
......@@ -18,14 +18,7 @@ package MixtureLiquid "Medium model of a mixture of liquids based on Modelica li
Density(start=10, nominal=10),
AbsolutePressure(start=10e5, nominal=10e5),
Temperature(min=273.15, max=573.15, start=300, nominal=300));
//substanceNames=data[:].name,
// fixedX=true,
// reducedX = false,
// singleState=false,
// nXi = nS-nR-1,
//
// SpecificEnthalpy(start=if Functions.referenceChoice==ReferenceEnthalpy.ZeroAt25C then 3e5 else
// if Functions.referenceChoice==ReferenceEnthalpy.UserDefined then Functions.h_offset else 0, nominal=1.0e5),
redeclare record extends ThermodynamicState "Thermodynamic state variables"
Real[nF] Xfull;
......@@ -40,16 +33,10 @@ package MixtureLiquid "Medium model of a mixture of liquids based on Modelica li
constant Integer nL = userInterface.nL "Number of liquid species (solutes+solvent) from user interface";
constant Integer nLi = nL-1 "Number of liquid species (solutes+solvent) -1";
constant SI.MassFraction[nL] reference_Xfull= userInterface.refXfull "Reference mass fraction vector of all species in the system";
//max(1,nF-1);
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= 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 MolarMass MH2O = IF97.MH2O "Molar mass of solvent";
......
......@@ -36,33 +36,16 @@ algorithm
z :=tau./n;
diagH[1:ns] :=tau ./ n[1:ns] .^ 2;//z[1:ns]./n[1:ns];//
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;
diagH[1:ns] :=tau ./ n[1:ns] .^ 2;
diagH[ns+1:ns+nL] :=(1 .- Yl .+ z[ns+1:ns+nL]) ./ n[ns + 1:ns + nL];
H :=diagonal(diagH);//fill(diagH,ns+nL);//
H :=diagonal(diagH);
H_id :=transpose(P_to_id)*H;
// 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;//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;
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);
J[nR+1:nF,:] :=transpose(lambda_id);
end SLE_Tp;
within ElectrolyteMedia.Media.SolidLiquidPhase.Common.MixtureLiquid.Initialization;
package calc_J
end calc_J;
......@@ -9,7 +9,6 @@ function SLE_Tp
protected
SI.MolarEnergy gi[ns+nL];
SI.MolarEnergy gi_id[ns+nL];
SI.MolarEnergy gr[nR];
Real[nR] logK;
Real[nL] gamma_i;
......@@ -22,24 +21,19 @@ protected
SI.MoleFraction[ns] Ys;
SI.MassFraction[ns] Xs;
Real tau = 1e-37;
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;
SI.MassFraction[ns+nL] z;
Real[ns+nL] logreacBase;
Real[ns+nL] logreacBase_id;
Real[ns+nL] x_orig;
Real[ns+nL] x_orig "moles in original order";
SI.Mass[ns+nL] mass;
algorithm
//assert(sum(x[1:ns+nL]) > 0, "x is zero in SLE_Tp");
assert(sum(x[1:ns+nL]) > 0, "x is zero in SLE_Tp");
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;
// z :={if x_orig[i] > 0 then tau/max(1e-20*tau, x_orig[i]) else 0 for i in 1:ns + nL};
z :=tau ./ x_orig;
z :=tau ./ x_orig;
Y :=x_orig[1:ns+nL]/sum(x_orig[1:ns+nL]);
X :=Functions.calc_Xfull(Y);
......@@ -55,24 +49,16 @@ algorithm
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*gi;
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]);
logreacBase_id :=P_to_id*logreacBase;
// x_id :=P_to_id*x;
//isopotential
f[1:nR] :=nu*logreacBase - logK;
// f[1:nR] :=nu_id*logreacBase_id - logK;
f[1:nR] :=nu*logreacBase - 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_id)*x - Xred;//transpose(lambda_mass)*mass - Xred;
f[nR+1:nF] :=transpose(lambda_id)*x - Xred;
end SLE_Tp;
within ElectrolyteMedia.Media.SolidLiquidPhase.Common.MixtureLiquid.Initialization;
package calc_f
end calc_f;
......@@ -10,6 +10,4 @@ package Initialization "Package containing all functions regarding initializatio
// constant DataRecordL[nL-1] datal;
// constant Media.Common.Types.LiquidModel LiquidModel;
end Initialization;
......@@ -31,7 +31,6 @@ protected
Integer index;
Real[nX] Xs;
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
......@@ -43,12 +42,12 @@ algorithm
Xfull :=cat(1,zeros(nF - nX),Xs);
Xfull :=P_to_orig*Xfull;
Xred_ :=transpose(lambda_mass)*Xfull;
Xred_ :=transpose(lambda_mass_id_orig)*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_id_orig)*Xinit);
minit :=Xinit*scale;
ninit :=minit./MMX;
for i in 1:nF loop
......
......@@ -21,7 +21,6 @@ package MixtureLiquid
redeclare record extends ThermodynamicState "Thermodynamic state variables"
Real[nF] Xfull;
end ThermodynamicState;
constant SI.Pressure prefg = 1e5 "Reference pressure for gas phase";
constant DataRecordL[nL-1] datal=userInterface.datal;
constant LiquidInteractionDataRecord interactionL=userInterface.interactionL;
......@@ -41,15 +40,15 @@ package MixtureLiquid
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,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 UserInterface
userInterface "User interface"
annotation (Placement(transformation(extent={{-10,-8},{10,12}})));
constant MolarMass MH2O = IF97.MH2O "Molar mass of solvent";
constant SpecificHeatCapacity RH2O = IF97.RH2O;
......@@ -156,9 +155,8 @@ package MixtureLiquid
logK = -gr / Modelica.Constants.R / T;
K = exp(logK);
log10K = log10(K);
logreacBase = cat(1,-z[1:ns],loga);//-z[1+ns:ns+nL]);
// nu*logreacBase = logK;
nu*logreacBase./logK = ones(nR);
logreacBase = cat(1,-z[1:ns],loga);
nu*logreacBase./logK = ones(nR);
a = exp(loga);
pH = -loga[Hindex-ns]*log10(exp(1));
Phil = Functions.calc_Phil(T,p,Xfull);
......@@ -267,7 +265,6 @@ protected
annotation(Inline = true, smoothOrder = 3);
end density;
redeclare function extends molarMass "Return molar mass of mixture"
algorithm
MM := Functions.calc_MM(state.Xfull);
......@@ -354,16 +351,4 @@ protected
eta := 1;//Common.Functions.IF97_R1_Tp.calc_eta(state.p,state.T);
annotation (smoothOrder=2);
end dynamicViscosity;
end MixtureLiquid;
......@@ -17,6 +17,8 @@ nu_id
P_to_orig
P_to_id
lambda_id
lambda_mass_id
lambda_mass_id_orig
MMX_id
userInterface
MH2O
......
......@@ -233,7 +233,7 @@ package SolidLiquidPhase
SI.MassFlowRate[ns+nL,N] m_out;
initial equation
for n in 1:N loop
mred_n[:,n] = transpose(lambda_mass)*medium_n[n].Xfullstart*d_start*B*H*deltaL;
mred_n[:,n] = transpose(lambda_mass)*Medium.userInterface.refXfull*d_start*B*H*deltaL;
end for;
equation
......@@ -292,8 +292,8 @@ package SolidLiquidPhase
lfrac_in=1)
annotation (Placement(transformation(extent={{-10,-10},{10,10}})));
annotation (experiment(
StopTime=60000,
Tolerance=1e-06,
__Dymola_Algorithm="Lsodar"));
StopTime=60000,
Tolerance=1e-06,
__Dymola_Algorithm="Lsodar"));
end ReactiveTransportSimulation;
end SolidLiquidPhase;
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment