I am working with OMEdit (v1.22.3) and I tried to replicate something similar to the model shown in this journal paper:
“Influence of design on performance of a latent heat storage system for a
direct steam generation power plant” (Website: https://www.sciencedirect.com/science/article/abs/pii/S030626191501332X)
It is about the discharge of a storage where there is a phase change material using subcooled water.
I was working with some functions to calculate water/steam properties, but I think there is something wrong. When I run the model “Module”, it calculates fair enough the solution until the first phase change happens. Later, it seems that the solution does not converge for any reason I don’t know. I have no clue what the problem can be, so If someone can help me, I will be very grateful. I share the code here.
package TFM_V0
package h_conv
function h_biphasic
input Modelica.Units.SI.Pressure P "Pressure of water";
input Modelica.Units.SI.QualityFactor x_HTF "Quality of water";
input Modelica.Units.SI.Temperature T_s "Surface temperature";
input Modelica.Units.SI.MassFlowRate m_flow "Mass flow rate";
input Modelica.Units.SI.Length D "Diameter of the pipe";
output Modelica.Units.SI.CoefficientOfHeatTransfer h_htf "Single-phase convection heat transfer coefficient";
protected
constant Real pi = 3.14159265359;
constant Real g(unit = "m/s2") = 9.81;
Modelica.Units.SI.SpecificEnthalpy h_l(start = 1e-3, fixed = false);
Modelica.Units.SI.SpecificEnthalpy Deltah_lv(start = 1e-3, fixed = false);
Modelica.Units.SI.Density rho_l(start = 1e-3, fixed = false);
Modelica.Units.SI.Density rho_v(start = 1e-3, fixed = false);
Modelica.Units.SI.Temperature T_sat(start = 1e-3, fixed = false);
Modelica.Units.SI.SurfaceTension sigma(start = 1e-3, fixed = false);
Modelica.Units.SI.DynamicViscosity mu_l(start = 1e-3, fixed = false);
Modelica.Units.SI.DynamicViscosity mu_v(start = 1e-3, fixed = false);
Modelica.Units.SI.ThermalConductivity k_l(start = 1e-3, fixed = false);
Modelica.Units.SI.SpecificHeatCapacity cp_l(start = 1e-3, fixed = false);
Modelica.Units.SI.Pressure P_sat_T_s(start = 1e-3, fixed = false);
//Dimensionless numbers
Real inv_X_tt;
Real Re_l;
Real Re_FB;
Real Pr_l;
Real F;
Real S;
Real x;
//Real C;
Modelica.Units.SI.CoefficientOfHeatTransfer h_htf_mac;
Modelica.Units.SI.CoefficientOfHeatTransfer h_htf_mic;
//Package properties of fluid: Water
package Fluid = Modelica.Media.Water.IF97_Utilities;
algorithm
//Correction quality factor
if x_HTF == 1 then
x := 0.99;
elseif x_HTF == 0 then
x := 0.01;
else
x := x_HTF;
end if;
//Thermo-physical properties
h_l := Fluid.hl_p(P);
Deltah_lv := Fluid.hv_p(P) - h_l;
rho_l := Fluid.rhol_p(P);
rho_v := Fluid.rhov_p(P);
T_sat := Fluid.BaseIF97.Basic.tsat(P);
sigma := Fluid.surfaceTension(T_sat);
mu_l := Fluid.dynamicViscosity(rho_l, T_sat, P);
mu_v := Fluid.dynamicViscosity(rho_v, T_sat, P);
k_l := Fluid.thermalConductivity(rho_l, T_sat, P);
cp_l := Fluid.cp_ph(P, h_l, 0, 0);
P_sat_T_s := Fluid.BaseIF97.Basic.psat(T_s);
//Dimensionless numbers
Pr_l := cp_l * mu_l / k_l;
Re_l := D * (1 - x) * (4 * m_flow / (pi * D ^ 2)) / mu_l;
inv_X_tt := (x / (1 - x)) ^ 0.9 * (rho_l / rho_v) ^ 0.5 * (mu_v / mu_l) ^ 0.1;
if inv_X_tt <= 0.1 then
F := 1;
else
F := 2.35 * (inv_X_tt + 0.213) ^ 0.736;
end if;
Re_FB := Re_l * F ^ 1.25 * 1e-4;
if Re_FB < 32.5 then
S := 1 / (1 + 0.12 * Re_FB ^ 1.14);
else
if Re_FB > 70 then
S := 0.1;
else
S := 1 / (1 + 0.42 * Re_FB ^ 0.78);
end if;
end if;
//Convective heat transfer coefficient
h_htf_mac := 0.023 * Re_l ^ 0.8 * Pr_l ^ 0.4 * k_l / D * F;
h_htf_mic := 0.00122 * (k_l ^ 0.79 * cp_l ^ 0.45 * rho_l ^ 0.49 * g ^ 0.25 / (sigma ^ 0.5 * mu_l ^ 0.29 * Deltah_lv ^ 0.24 * rho_v ^ 0.24)) * abs(T_s - T_sat) ^ 0.24 * abs(P_sat_T_s - P) ^ 0.75 * S;
h_htf := h_htf_mic + h_htf_mac;
end h_biphasic;
function h_single
input Modelica.Units.SI.Pressure P "Pressure of water";
input Modelica.Units.SI.Temperature T "Temperature of water";
input Modelica.Units.SI.MassFlowRate m_flow "Mass flow rate";
input Modelica.Units.SI.Length D "Diameter of the pipe";
output Modelica.Units.SI.CoefficientOfHeatTransfer h_htf "Single-phase convection heat transfer coefficient";
protected
constant Real pi = 3.14159265359;
//Thermo-physical properties
Modelica.Units.SI.Density rho(start = 1e-3, fixed = false);
Modelica.Units.SI.SpecificHeatCapacity cp(start = 1e-3, fixed = false);
Modelica.Units.SI.SpecificEnthalpy h(start = 1e-3, fixed = false);
Modelica.Units.SI.ThermalConductivity k(start = 1e-3, fixed = false);
Modelica.Units.SI.DynamicViscosity mu(start = 1e-3, fixed = false);
//Dimensionless numbers
Real Re;
Real Pr;
Real f;
Real Nu;
//Package properties of fluid: Water
package Fluid = Modelica.Media.Water.IF97_Utilities;
algorithm
//Thermal properties of single phase
h := Fluid.h_pT(P, T, 0);
rho := Fluid.rho_ph(P, h, 0, 0);
cp := Fluid.cp_ph(P, h, 0, 0);
k := Fluid.thermalConductivity(rho, T, P);
mu := Fluid.dynamicViscosity(rho, T, P);
//Dimensionless numbers
Pr := cp * mu / k;
Re := 4 * m_flow / (pi * D * mu);
f := (1.82 * log10(Re) - 1.64) ^ (-2);
Nu := f / 8 * (Re - 1000) * Pr / (1 + 12.7 * (f / 8) ^ 0.5 * (Pr ^ (2 / 3) - 1));
//Convective heat transfer coefficient
h_htf := Nu * k / D;
end h_single;
end h_conv;
package Two_phase_props
function h_T
input Modelica.Units.SI.Temperature T;
input Modelica.Units.SI.SpecificHeatCapacity cp_low;
input Modelica.Units.SI.SpecificHeatCapacity cp_high;
input Modelica.Units.SI.SpecificEnthalpy h_low;
input Modelica.Units.SI.SpecificEnthalpy h_high;
input Modelica.Units.SI.Temperature T_phasechange;
output Modelica.Units.SI.SpecificEnthalpy h;
algorithm
if T > T_phasechange then
h := h_high + cp_high * (T - T_phasechange);
else
h := h_low + cp_low * (T - T_phasechange);
end if;
end h_T;
function T_h
input Modelica.Units.SI.SpecificEnthalpy h;
input Modelica.Units.SI.SpecificHeatCapacity cp_low;
input Modelica.Units.SI.SpecificHeatCapacity cp_high;
input Modelica.Units.SI.SpecificEnthalpy h_low;
input Modelica.Units.SI.SpecificEnthalpy h_high;
input Modelica.Units.SI.Temperature T_phasechange;
output Modelica.Units.SI.Temperature T;
algorithm
if h >= h_low and h <= h_high then
T := T_phasechange;
else
if h > h_high then
T := T_phasechange + (h - h_high) / cp_high;
else
T := T_phasechange + (h - h_low)/ cp_low;
end if;
end if;
end T_h;
function x_h
input Modelica.Units.SI.SpecificEnthalpy h;
input Modelica.Units.SI.SpecificEnthalpy h_low;
input Modelica.Units.SI.SpecificEnthalpy h_high;
output Modelica.Units.SI.QualityFactor x "Quality of HTF";
algorithm
if h <= h_high and h >= h_low then
x := (h - h_low) / (h_high - h_low);
else
if h >= h_low then
x := 1;
else
x := 0;
end if;
end if;
end x_h;
end Two_phase_props;
model Module
//Constants
constant Real pi = 3.14159265359;
constant Real g (unit = "m/s2") = 9.81;
//Design of module (Geometry)
parameter Modelica.Units.SI.Length Length = 100 "Length of the module";
parameter Modelica.Units.SI.Length d_int = 0.04094 "Inner pipe diameter";
parameter Modelica.Units.SI.Length D_ext = 0.0483 "External pipe diameter";
parameter Modelica.Units.SI.Length D_mod = 0.5 "Diameter of the module";
//Initial conditions and constant operation
parameter Modelica.Units.SI.MassFlowRate m_flow_nom = 1 / 3.6 "HTF mass flow";
parameter Modelica.Units.SI.Pressure P_nom = 1000000 "HTF Nominal pressure";
parameter Real T_discharge(unit = "K", displayUnit = "degC") = 373.15 "Low dischrage temperature";
parameter Modelica.Units.SI.QualityFactor x_discharge = 0 "Low dischrage quality";
//Pipe
parameter Modelica.Units.SI.ThermalConductivity k_pipe = 21 "Thermal conductivity of pipe";
//Module properties
parameter Real T_ini(unit = "K", displayUnit = "degC") = 573.15 "PCM initial temperature";
parameter Modelica.Units.SI.Density rho_s_PCM = 2500 "Solid-phase density";
parameter Modelica.Units.SI.Density rho_l_PCM = 2300 "Liquid-phase density";
parameter Modelica.Units.SI.SpecificHeatCapacity cp_s_PCM = 900 "Solid-phase specific heat capacity";
parameter Modelica.Units.SI.SpecificHeatCapacity cp_l_PCM = 1600 "Liquid-phase specific heat capacity";
parameter Modelica.Units.SI.SpecificEnthalpy L_PCM = 120000 "PCM latent heat of fusion";
parameter Real T_melt(unit = "K", displayUnit = "degC") = 478.15 "Melting temperature";
parameter Modelica.Units.SI.ThermalConductivity k_PCM = 2.2 "Thermal conductivity of PCM";
parameter Modelica.Units.SI.Mass m_PCM = rho_l_PCM * pi / 4 * (D_mod ^ 2 - D_ext ^ 2) "Total mass of PCM in the module";
//HTF properties
parameter Real T_sat_HTF (unit = "K", displayUnit = "degC", start = 1e-3, fixed = false) "HTF saturation tempearture";
//Discretization
parameter Integer nodes = 20 "Number of nodes";
parameter Modelica.Units.SI.Volume v_HTF_node = pi / 4 * d_int ^ 2 * dL;
//Others
parameter Modelica.Units.SI.Length dL = Length / nodes "Slice thickness";
//Variables
//A: HTF
Modelica.Units.SI.SpecificEnthalpy h_HTF[nodes] (each start = 3.05e6) "HTF enthalpy in slices";
Real T_HTF[nodes] (each start = 573.15, each unit = "K", each displayUnit = "degC") "Temperature of HTF slices";
Modelica.Units.SI.QualityFactor x_HTF[nodes](each min = 0, each max = 1) "Quality of HTF slices";
Modelica.Units.SI.Density rho_HTF[nodes] "Density of HTF at each slice";
Modelica.Units.SI.SpecificEnthalpy h_l "Saturated liquid enthalpy at Pnom";
Modelica.Units.SI.SpecificEnthalpy h_v "Saturated vapor enthalpy at Pnom";
Modelica.Units.SI.SpecificEnthalpy h_discharge "Enthalpy at discharge";
//B: PCM
Modelica.Units.SI.SpecificEnthalpy h_PCM[nodes](each start = 272e3) "PCM enthalpy in slices";
Real T_PCM[nodes] (each start = 573.15, each unit = "K", each displayUnit = "degC") "Temperature of PCM slices";
Modelica.Units.SI.QualityFactor x_PCM[nodes](each min = 0, each max = 1) "Melting fraction of PCM slices";
//C: Module
Real T_si[nodes] (each start = 573.15, each unit = "K", each displayUnit = "degC") "Internal wall temperature of slices";
Real T_so[nodes] (each start = 573.15, each unit = "K", each displayUnit = "degC") "External wall temperature of slices";
Modelica.Units.SI.CoefficientOfHeatTransfer h_conv[nodes] "Convective heat transfer coefficient in slices";
Real UA[nodes](each unit = "W/K") "UA product in slices";
Modelica.Units.SI.HeatFlowRate Q_dot[nodes] "Heat flow in slices";
Modelica.Units.SI.Length D_PCM[nodes] "Diameter of the solid PCM in slices";
Real T_outlet (start = 573.15, unit = "K", displayUnit = "degC") "HTF outlet temperature";
//Packages: Properties of water
replaceable package HTF = Modelica.Media.Water.IF97_Utilities;
//Functions
replaceable function h_conv_SP = TFM_V0.h_conv.h_single;
replaceable function h_conv_TP = TFM_V0.h_conv.h_biphasic;
replaceable function h_T = TFM_V0.Two_phase_props.h_T;
replaceable function T_h = TFM_V0.Two_phase_props.T_h;
replaceable function x_h = TFM_V0.Two_phase_props.x_h;
initial equation
for i in 1:nodes loop
//PCM
T_PCM[i] = T_ini;
h_PCM[i] = h_T(T_ini,cp_s_PCM,cp_l_PCM,0,L_PCM,T_melt);
x_PCM[i] = x_h(h_PCM[i],0,L_PCM);
end for;
//HTF: Saturation properties
T_sat_HTF = HTF.BaseIF97.Basic.tsat(P_nom);
//HTF: initial properties
x_discharge = x_h(h_discharge,h_l,h_v);
equation
//HTF: Saturation properties
h_l = HTF.hl_p(P_nom);
h_v = HTF.hv_p(P_nom);
//HTF: initial properties
h_discharge = HTF.h_pT(P_nom,T_discharge,1);
for i in 1:nodes loop
//HTF conditions
x_HTF[i] = noEvent(if (h_HTF[i] < h_v and h_HTF[i] > h_l) then (h_HTF[i] - h_l) / (h_v - h_l)
else (if h_HTF[i] >= h_v then 1 else 0));
T_HTF[i] = smooth(0, HTF.T_ph(P_nom,h_HTF[i],0,0));
rho_HTF[i] = smooth(0, HTF.rho_ph(P_nom,h_HTF[i],0,0));
//Slice model in HTF
if i == 1 then
Q_dot[1] = rho_HTF[1] * v_HTF_node * der(h_HTF[1]) + m_flow_nom * (h_HTF[1] - h_discharge);
else
Q_dot[i] = rho_HTF[i] * v_HTF_node * der(h_HTF[i]) + m_flow_nom * (h_HTF[i] - h_HTF[i-1]);
end if;
//PCM conditions: Heat flow is "negative" because the PCM "loss" heat
x_PCM[i] = noEvent(if (h_PCM[i] < L_PCM and h_PCM[i] > 0) then h_PCM[i] / L_PCM
else (if h_PCM[i] >= L_PCM then 1 else 0));
T_PCM[i] = noEvent(T_h(h_PCM[i],cp_s_PCM,cp_l_PCM,0,L_PCM,T_melt));
der(h_PCM[i]) = -Q_dot[i] / (m_PCM / nodes);
//Convection heat transfer
if x_HTF[i] > 0 and x_HTF[i] < 1 then
h_conv[i] = h_conv_TP(P_nom,x_HTF[i],T_si[i],m_flow_nom,d_int);
else
h_conv[i] = h_conv_SP(P_nom,T_HTF[i],m_flow_nom,d_int);
end if;
//Heat transfer through pipe: Convection heat flow = Conduction heat flow in the pipe = Heat flow from PCM to HTF
h_conv[i] * pi * d_int * dL * (T_si[i] - T_HTF[i]) = 2 * pi * k_pipe * dL * (T_so[i] - T_si[i]) / log(D_ext / d_int);
Q_dot[i] = 2 * pi * k_pipe * dL * (T_so[i] - T_si[i]) / log(D_ext / d_int);
//Space of solid-phase PCM in a slice
D_PCM[i] ^ 2 = (1 - x_PCM[i]) * (D_mod ^ 2 - D_ext ^ 2) + D_ext ^ 2;
//Slice model
Q_dot[i] = UA[i] * (T_PCM[i] - T_HTF[i]);
if x_PCM[i] > 0 and x_PCM[i] < 1 then
UA[i] = 1 / (1 / (h_conv[i] * pi * d_int * dL) + log(D_ext / d_int) / (2 * pi * k_pipe * dL) +
log(D_PCM[i] / D_ext) / (2 * pi * k_PCM * dL));
else
UA[i] = 1 / (1 / (h_conv[i] * pi * d_int * dL) + log(D_ext / d_int) / (2 * pi * k_pipe * dL));
end if;
end for;
T_outlet = T_HTF[nodes];
end Module;
end TFM_V0;
I tried noEvent at every if-cycle in the model, also smooth() (there are still some that has worked “fair enough” as I said before).