Skip to content

Commit

Permalink
reimplementation of discretized hex and refactor of volumes (#31)
Browse files Browse the repository at this point in the history
* mirnor bugfixes and refactor of volume models
* add volume models with multiple outlets and phase separator with two outlets
* reimplementation of Discretized HEX's with conductionElements and addition of cross-flow variants
* refactor of HEX tests
* minor bugfixes
  • Loading branch information
mimeissner authored Nov 17, 2021
1 parent c9839ea commit f2eb7b0
Show file tree
Hide file tree
Showing 59 changed files with 3,891 additions and 1,437 deletions.
3 changes: 1 addition & 2 deletions ThermofluidStream/Boundaries/Internal/PartialVolume.mo
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,6 @@ protected
Medium.MassFraction Xi_in[Medium.nXi] = if noEvent(m_flow_in >= 0) then Medium.massFraction(state_in) else medium.Xi;

Medium.ThermodynamicState state_out;
SI.Pressure p_out = Medium.pressure(state_out);
// fix potential instabilities by setting the incoming enthalpy and mass fraction inlet ones,
// effectiveley removing the mass-flow related parts of the differential equations for U and MXi
SI.SpecificEnthalpy h_out = if noEvent(-m_flow_out >= 0) then Medium.specificEnthalpy(state_out) else medium.h;
Expand Down Expand Up @@ -112,7 +111,7 @@ equation
der(m_flow_in)*L = r_in - r - r_damping;
der(m_flow_out)*L = r_out - r_damping;

r + p_in = p_out;
r + p_in = medium.p;

der(M) = m_flow_in + m_flow_out;
der(U_med) = W_v + Q_flow + h_in*m_flow_in + h_out*m_flow_out;
Expand Down
193 changes: 193 additions & 0 deletions ThermofluidStream/Boundaries/Internal/PartialVolumeM.mo
Original file line number Diff line number Diff line change
@@ -0,0 +1,193 @@
within ThermofluidStream.Boundaries.Internal;
partial model PartialVolumeM "Partial parent class for Volumes with one inlet and M outlet"
replaceable package Medium = Media.myMedia.Interfaces.PartialMedium "Medium model" annotation (
choicesAllMatching=true, Documentation(info="<html>
<p><span style=\"font-family: Courier New;\">Medium package used in the Volume. Make sure it is the same as the inlets and outlets the volume is connected to.</span></p>
</html>"));

parameter Integer M_outlets = 1 "Number if outlets";
parameter Boolean useHeatport = false "If true heatport is added";
parameter SI.Area A = 1 "Contact area of volume with medium"
annotation(Dialog(enable=useHeatport));
parameter SI.CoefficientOfHeatTransfer U = 200 "Heat transfer coefficient to medium"
annotation(Dialog(enable=useHeatport));
parameter Boolean initialize_pressure = true "If true: initialize Pressure"
annotation(Dialog(tab= "Initialization"));
parameter SI.Pressure p_start = Medium.p_default "Initial Pressure"
annotation(Dialog(tab= "Initialization", enable=initialize_pressure));
parameter Boolean initialize_energy = true "Initialize specific inner energy with temperature or specific enthalpy condition"
annotation(Dialog(tab= "Initialization"));
parameter SI.Temperature T_start = Medium.T_default "Initial Temperature"
annotation(Dialog(tab= "Initialization", enable=initialize_energy and (not use_hstart)));
parameter Boolean use_hstart = false "True: spedific enthalpy contition instead of Temperature"
annotation(Dialog(tab= "Initialization", enable=initialize_energy));
parameter SI.SpecificEnthalpy h_start = Medium.T_default "Initial specific enthalpy"
annotation(Dialog(tab= "Initialization", enable=initialize_energy and use_hstart));
parameter Boolean initialize_Xi = true "If true: initialize mass fractions"
annotation(Dialog(tab= "Initialization"));
parameter Medium.MassFraction Xi_0[Medium.nXi] = Medium.X_default[1:Medium.nXi] "Initial mass fraction"
annotation(Dialog(tab= "Initialization", enable=initialize_Xi));
parameter Utilities.Units.Inertance L = dropOfCommons.L "Inertance at inlet and outlet"
annotation (Dialog(tab="Advanced"));
parameter Real k_volume_damping(unit="1") = dropOfCommons.k_volume_damping "Damping factor multiplicator"
annotation(Dialog(tab="Advanced", group="Damping"));
parameter SI.MassFlowRate m_flow_assert(max=0) = -dropOfCommons.m_flow_reg "Assertion threshold for negative massflows"
annotation(Dialog(tab="Advanced"));
parameter Boolean usePreferredMediumStates=false "Use medium states instead of the ones differentiated in this component"
annotation(Dialog(tab="Advanced"));

Interfaces.Inlet inlet(redeclare package Medium=Medium)
annotation (Placement(transformation(extent={{-120,-20},{-80,20}})));
Interfaces.Outlet outlet[M_outlets](redeclare package Medium=Medium)
annotation (Placement(transformation(extent={{80,-20},{120,20}})));
Modelica.Thermal.HeatTransfer.Interfaces.HeatPort_a heatPort(Q_flow=Q_flow, T=T_heatPort) if useHeatport
annotation (Placement(transformation(extent={{-10,-90},{10,-70}})));

Medium.BaseProperties medium(preferredMediumStates=usePreferredMediumStates);

SI.Volume V;

//setting the state is to prohibit dynamic state selection e.g. in VolumesDirectCoupling
SI.Mass M(stateSelect=if usePreferredMediumStates then StateSelect.default else StateSelect.always) = V*medium.d;
SI.Mass MXi[Medium.nXi](each stateSelect=if usePreferredMediumStates then StateSelect.default else StateSelect.always) = M*medium.Xi;
SI.Energy U_med(stateSelect=if usePreferredMediumStates then StateSelect.default else StateSelect.always) = M*medium.u;

SI.HeatFlowRate Q_flow;
SI.Power W_v;

protected
outer DropOfCommons dropOfCommons;

SI.Pressure p_in = Medium.pressure(inlet.state);
// fix potential instabilities by setting the outgoing enthalpy and mass fraction to the medium state
SI.SpecificEnthalpy h_in = if noEvent(m_flow_in) >= 0 then Medium.specificEnthalpy(inlet.state) else medium.h;
Medium.MassFraction Xi_in[Medium.nXi] = if noEvent(m_flow_in) >= 0 then Medium.massFraction(inlet.state) else medium.Xi;

Medium.ThermodynamicState state_out[M_outlets];
SI.SpecificEnthalpy h_out[M_outlets];
Medium.MassFraction Xi_out[Medium.nXi,M_outlets];

Real d(unit="1/(m.s)") = k_volume_damping*sqrt(abs(2*L/(V*max(density_derp_h, 1e-10)))) "Friction factor for coupled boundaries";
SI.DerDensityByPressure density_derp_h "Partial derivative of density by pressure";
SI.Pressure r_damping = d*der(M);

SI.Pressure r;

SI.Temperature T_heatPort;

SI.MassFlowRate m_flow_in = inlet.m_flow;
SI.MassFlowRate m_flow_out[M_outlets] = outlet.m_flow;

initial equation
if initialize_pressure then
medium.p=p_start;
end if;

if initialize_energy then
if use_hstart then
medium.h = h_start;
else
medium.T=T_start;
end if;
end if;

if initialize_Xi then
medium.Xi = Xi_0;
end if;

equation
assert(m_flow_in > m_flow_assert, "Negative massflow at Volume inlet", dropOfCommons.assertionLevel);
for i in 1:M_outlets loop
assert(-m_flow_out[i] > m_flow_assert, "Positive massflow at Volume outlet", dropOfCommons.assertionLevel);
end for;
assert(M > 0, "Volumes might not become empty");

der(inlet.m_flow)*L = inlet.r - r - r_damping;
der(outlet.m_flow)*L = outlet.r - r_damping*ones(M_outlets);

r + p_in = medium.p;

for i in 1:M_outlets loop
// fix potential instabilities by setting the outgoing enthalpy and mass fraction to the medium state

h_out[i] = if noEvent(-m_flow_out[i]) >= 0 then Medium.specificEnthalpy(state_out[i]) else medium.h;
Xi_out[:,i] = if noEvent(-m_flow_out[i] >= 0) then Medium.massFraction(state_out[i]) else medium.Xi;
end for;

der(M) = inlet.m_flow + sum(outlet.m_flow);
der(U_med) = W_v + Q_flow + inlet.m_flow*h_in + sum(outlet.m_flow.*h_out);
der(MXi) = Xi_in*inlet.m_flow + Xi_out*outlet.m_flow;

Q_flow = U*A*(T_heatPort - medium.T);

outlet.state = state_out;

if not useHeatport then
T_heatPort = medium.T;
end if;

annotation (Icon(coordinateSystem(preserveAspectRatio=false), graphics={
Ellipse(
extent={{-56,76},{64,16}},
lineColor={28,108,200},
lineThickness=0.5,
fillColor={215,215,215},
fillPattern=FillPattern.Solid,
pattern=LinePattern.None),
Rectangle(
extent={{-56,46},{64,-56}},
lineColor={28,108,200},
lineThickness=0.5,
fillColor={215,215,215},
fillPattern=FillPattern.Solid,
pattern=LinePattern.None),
Ellipse(
extent={{-56,-28},{64,-88}},
lineColor={28,108,200},
lineThickness=0.5,
fillColor={215,215,215},
fillPattern=FillPattern.Solid,
pattern=LinePattern.None),
Line(
points={{-100,0},{100,0}},
color={28,108,200},
thickness=0.5),
Ellipse(
extent={{-60,-20},{60,-80}},
lineColor={28,108,200},
lineThickness=0.5,
fillColor={170,213,255},
fillPattern=FillPattern.Solid),
Rectangle(
extent={{-60,50},{60,-50}},
lineColor={28,108,200},
lineThickness=0.5,
fillColor={170,213,255},
fillPattern=FillPattern.Solid,
pattern=LinePattern.None),
Ellipse(
extent={{-60,80},{60,20}},
lineColor={28,108,200},
lineThickness=0.5,
fillColor={170,213,255},
fillPattern=FillPattern.Solid),
Line(
points={{-60,50},{-60,-52}},
color={28,108,200},
thickness=0.5),
Line(
points={{60,50},{60,-52}},
color={28,108,200},
thickness=0.5), Text(
extent={{68,48},{94,6}},
lineColor={116,116,116},
textString="%M_out")}),
Diagram(coordinateSystem(preserveAspectRatio=false)),
Documentation(info="<html>
<p>This is the partial parent class for unidirectional volumes with N inlets. It is partial and missing the number if inputs, as well as a equations for its volume or medium pressure and one for the volume change work performed.</p>
<p>Conceptually&nbsp;a&nbsp;Volume&nbsp;is&nbsp;a&nbsp;Sink&nbsp;and&nbsp;a&nbsp;Source.&nbsp;It&nbsp;therefore&nbsp;defines&nbsp;the&nbsp;Level&nbsp;of&nbsp;inertial&nbsp;pressure&nbsp;r&nbsp;in&nbsp;a&nbsp;closed&nbsp;loop&nbsp;and&nbsp;acts&nbsp;as&nbsp;a&nbsp;Loop&nbsp;breaker.</p>
<p>Volumes implement a damping term on the change of the stored mass to dampen out fast, otherwise undamped oscillations that appear when connecting volumes directly to other volumes or other boundaries (source, sink, boundary_fore, boundary_rear). With the damping term these oscillations will be still very fast, but dampend out, so a stiff solver might be able to handle them well. Damping is enabled by default and can be disabled by setting Advanced.k_volume_damping=0. </p>
<p>For to stability reasons, mass-flows in the wrong direction (fluid entering the outlet or exiting the inlet) is considered to have the enthalpy and mass-fractions of the medium in the volume. This results in a stable steady-state solution, since this method effectiveley removes the parts of the energy and mass-fraction differential equations, that are associated with mass-flows. </p>
<p>Per default the Volume has the two states energy and mass (U_med and M) and one state for each mass, as well as one state for each substance of the fluid (except the first one). These will be enforced to be states of the simulation, which can result in nonlinear systems of size one or two, but works very reliable. To get rid of these Systems the modeler can enable the flag &apos;usePreferredMediumStates&apos; in the &apos;Advanced&apos; tab. Then the volume uses the states prefered by the medium object, rather then the default ones, which can improve the nonlinear systems most of the time, but also might lead to larger nonlinear systems (e.g. in the Test &apos;VolumesDirectCoupling&apos;).</p>
</html>"));
end PartialVolumeM;
14 changes: 8 additions & 6 deletions ThermofluidStream/Boundaries/Internal/PartialVolumeN.mo
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ partial model PartialVolumeN "Partial parent class for Volumes with N inlets and
<p><span style=\"font-family: Courier New;\">Medium package used in the Volume. Make sure it is the same as the inlets and outlets the volume is connected to.</span></p>
</html>"));

parameter Integer N = 1 "Number if inputs";
parameter Integer N = 1 "Number if inlets";
parameter Boolean useHeatport = false "If true heatport is added";
parameter SI.Area A = 1 "Contact area of volume with medium"
annotation(Dialog(enable=useHeatport));
Expand Down Expand Up @@ -63,8 +63,10 @@ protected
SI.SpecificEnthalpy h_in[N];
Medium.MassFraction Xi_in[Medium.nXi,N];

SI.SpecificEnthalpy h_out = if noEvent(-m_flow_out) >= 0 then Medium.specificEnthalpy(medium.state) else medium.h;
Medium.MassFraction Xi_out[Medium.nXi] = if noEvent(-m_flow_out >= 0) then Medium.massFraction(medium.state) else medium.Xi;
Medium.ThermodynamicState state_out;
// fix potential instabilities by setting the outgoing enthalpy and mass fraction to the medium state
SI.SpecificEnthalpy h_out = if noEvent(-m_flow_out) >= 0 then Medium.specificEnthalpy(state_out) else medium.h;
Medium.MassFraction Xi_out[Medium.nXi] = if noEvent(-m_flow_out >= 0) then Medium.massFraction(state_out) else medium.Xi;

Real d(unit="1/(m.s)") = k_volume_damping*sqrt(abs(2*L/(V*max(density_derp_h, 1e-10)))) "Friction factor for coupled boundaries";
SI.DerDensityByPressure density_derp_h "Partial derivative of density by pressure";
Expand Down Expand Up @@ -113,12 +115,12 @@ equation
end for;

der(M) = sum(inlet.m_flow) + sum(outlet.m_flow);
der(U_med) = W_v + Q_flow + sum(inlet.m_flow.*h_in) + outlet.m_flow*medium.h;
der(MXi) = Xi_in*inlet.m_flow + medium.Xi*outlet.m_flow;
der(U_med) = W_v + Q_flow + sum(inlet.m_flow.*h_in) + outlet.m_flow*h_out;
der(MXi) = Xi_in*inlet.m_flow + Xi_out*outlet.m_flow;

Q_flow = U*A*(T_heatPort - medium.T);

outlet.state = medium.state;
outlet.state = state_out;

if not useHeatport then
T_heatPort = medium.T;
Expand Down
1 change: 1 addition & 0 deletions ThermofluidStream/Boundaries/Internal/package.order
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
PartialVolume
PartialVolumeN
PartialVolumeM
InitializationMethodsPhaseSeperator
99 changes: 99 additions & 0 deletions ThermofluidStream/Boundaries/PhaseSeparator2.mo
Original file line number Diff line number Diff line change
@@ -0,0 +1,99 @@
within ThermofluidStream.Boundaries;
model PhaseSeparator2 "Phase separator with two outlets"
extends Boundaries.Internal.PartialVolumeM(
redeclare replaceable package Medium =
Media.myMedia.Interfaces.PartialTwoPhaseMedium,
useHeatport=false,
final initialize_energy=false,
final T_start=0,
final h_start=0,
final use_hstart=false,
final M_outlets=2);

import Init = ThermofluidStream.Boundaries.Internal.InitializationMethodsPhaseSeperator;

parameter SI.Volume V_par(displayUnit="l")=0.01 "Volume of phase seperator";
parameter Real pipe1_low(unit="1", min=0, max=1) "Low end of pipe";
parameter Real pipe1_high(unit="1", min=0, max=1) "High end of pipe";
parameter Real pipe2_low(unit="1", min=0, max=1) "Low end of pipe";
parameter Real pipe2_high(unit="1", min=0, max=1) "High end of pipe";
parameter Boolean density_derp_h_from_media=false "EXPERIMENTAL: get density_derp_h from media model. The function is only implemented for some Media."
annotation(Dialog(tab="Advanced", group="Damping", enable=(k_volume_damping > 0)));
parameter SI.DerDensityByPressure density_derp_h_set = 1e-6 "Derivative of density by pressure upper bound; Approx. 1e-5 for air, 1e-7 for water"
annotation(Dialog(enable = ((k_volume_damping > 0) and not density_derp_h_from_media), tab="Advanced", group="Damping"));
parameter Init init_method = ThermofluidStream.Boundaries.Internal.InitializationMethodsPhaseSeperator.l "Initialization Method"
annotation(choicesAllMatching=true, Dialog(tab="Initialization"));
parameter SI.SpecificEnthalpy h_0 = Medium.h_default "Initial specific enthalpy"
annotation(Dialog(tab="Initialization", enable=(init_method==Init.h)));
parameter SI.Mass M_0 = 1 "Initial Mass"
annotation(Dialog(tab="Initialization", enable=(init_method==Init.M)));
parameter Real l_0(unit="1", min=0, max=1) = 0.5 "Initial liquid level"
annotation(Dialog(tab="Initialization", enable=(init_method==Init.l)));
parameter Real x_0(unit="kg/kg", min=0, max=1) = 0.5 "Initial vapor quality"
annotation(Dialog(tab="Initialization", enable=(init_method==Init.x)));

//Variables to calculate h_pipe
Real liquid_level(unit="1") "level of liquid line";
Real liquid_level_pipe1(unit="1") "level of liquid line in pipe1";
Real liquid_level_pipe2(unit="1") "level of liquid line in pipe2";

protected
Medium.MassFraction x = (medium.h-h_bubble)/(h_dew - h_bubble) "Calculated vapor quality of medium that can go below zero and above one";
SI.SpecificEnthalpy h_pipe1;
SI.SpecificEnthalpy h_pipe2;

SI.Density d_liq = Medium.bubbleDensity(Medium.setSat_p(medium.p)) "bubble density at saturation";
SI.Density d_gas = Medium.dewDensity(Medium.setSat_p(medium.p)) "dew density at saturation";

SI.SpecificEnthalpy h_bubble = Medium.bubbleEnthalpy(Medium.setSat_p(medium.p))-1 "Bubble Enthalpy of Medium";
SI.SpecificEnthalpy h_dew = Medium.dewEnthalpy(Medium.setSat_p(medium.p))+1 "Dew Enthalpy of Medium";

initial equation
assert(pipe1_high > pipe1_low, "Upper pipe end must be higher then lower end.", AssertionLevel.error);
assert(pipe2_high > pipe2_low, "Upper pipe end must be higher then lower end.", AssertionLevel.error);

if init_method == Init.h then
medium.h = h_0;
elseif init_method == Init.M then
x/d_gas+(1-x)/d_liq = V/M_0;
assert(x>=0 and x<=1, "Initialization by Mass might be inaccurate outside the two-phase region", AssertionLevel.warning);
elseif init_method == Init.l then
x = (d_gas*(1-l_0))/(d_liq*l_0+d_gas*(1-l_0));
elseif init_method == Init.x then
x = x_0;
end if;

equation
if density_derp_h_from_media then
density_derp_h = Medium.density_derp_h(medium.state);
else
density_derp_h = density_derp_h_set;
end if;

V = V_par;
W_v = 0;

liquid_level = max(0, min(1, M*(1-x)/d_liq/V));
liquid_level_pipe1 = max(0, min(1, (liquid_level-pipe1_low)/(pipe1_high-pipe1_low)));
liquid_level_pipe2 = max(0, min(1, (liquid_level-pipe2_low)/(pipe2_high-pipe2_low)));

h_pipe1 = smooth(1,
if x < 0 then medium.h
elseif x <= 1 then liquid_level_pipe1*h_bubble + (1-liquid_level_pipe1)*h_dew
else medium.h);

h_pipe2 = smooth(1,
if x < 0 then medium.h
elseif x <= 1 then liquid_level_pipe2*h_bubble + (1-liquid_level_pipe2)*h_dew
else medium.h);

state_out[1] = Medium.setState_phX(medium.p, h_pipe1, medium.Xi);
state_out[2] = Medium.setState_phX(medium.p, h_pipe2, medium.Xi);

annotation (Icon(coordinateSystem(preserveAspectRatio=false)),
Diagram(coordinateSystem(preserveAspectRatio=false)),
Documentation(info="<html>
<p>This Volume is the parent class for Accumulator and Receiver models that seperate the two phases and are able and output gas, liquid or two-phase medium, dependin on its liquid level and the height of the outlet. </p>
<p>Since there is no formula to compute density_derp_h for this volume, a upper bound has to be set in the parameter density_derp_h_set. Alternativeley the derivative can be taken from the media model for all the media that implement the corresponding forumla by setting density_derp_h_from_media=true (default:false).</p>
</html>"));
end PhaseSeparator2;
2 changes: 2 additions & 0 deletions ThermofluidStream/Boundaries/VolumeMix.mo
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,8 @@ equation
V = V_par;
W_v = 0;

state_out = medium.state;

annotation (Icon(coordinateSystem(preserveAspectRatio=false), graphics={
Text(
extent={{-60,8},{60,-52}},
Expand Down
1 change: 1 addition & 0 deletions ThermofluidStream/Boundaries/package.order
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@ Volume
VolumeFlex
VolumeMix
PhaseSeparator
PhaseSeparator2
Reservoir
CreateState
Tests
Expand Down
Loading

0 comments on commit f2eb7b0

Please sign in to comment.