Salt simulation elements

src/rtctools_channel_flow/modelica/Deltares/ChannelFlow/Salt/Elements/NodeSalty.mo

@@ -0,0 +1,20 @@
+within Deltares.ChannelFlow.Salt.Elements;
+
+model NodeSalty 
+  /*
+  This block is designed to be used together with the "salt_simulation_mixin" to calculate dispersive and advective transport
+   between salty reservoir elements, do not user in optimization.
+  */
+  extends Deltares.ChannelFlow.Internal.HQTwoPort;
+  extends QMForcing;
+
+equation
+
+  HQUp.Q + HQDown.Q + sum(QForcing) = 0;
+  HQUp.M[1] + HQDown.M[1] + sum(MForcing)  = 0; //Mass balance, can be used for one substance, like salt.
+
+  HQDown.H = HQUp.H;
+  HQUp.C = HQDown.C;
+
+  annotation(Icon(coordinateSystem(extent = {{-100, -100}, {100, 100}}, preserveAspectRatio = true, initialScale = 0.1, grid = {10, 10}), graphics = {Polygon(visible = true, fillColor = {255, 170, 0}, fillPattern = FillPattern.Solid, points = {{0, 50}, {-30, 40}, {30, -40}, {0, -50}, {-30, -40}, {30, 40}}), Polygon(visible = true, fillColor = {255, 0, 0}, fillPattern = FillPattern.Solid, points = {{-50, 0}, {-40, 30}, {-30, 40}, {30, -40}, {40, -30}, {50, 0}, {40, 30}, {30, 40}, {-30, -40}, {-40, -30}})}));
+end NodeSalty;

src/rtctools_channel_flow/modelica/Deltares/ChannelFlow/Salt/Elements/QMForcing.mo

@@ -0,0 +1,11 @@
+within Deltares.ChannelFlow.Salt.Elements;
+
+class QMForcing
+  import SI = Modelica.Units.SI;
+  parameter Integer n_QForcing(min = 0) = 0;
+
+  input SI.VolumeFlowRate QForcing[n_QForcing];
+  input SI.VolumeFlowRate MForcing[n_QForcing,1];
+
+  annotation(Icon(graphics = {Line(origin = {-40, 0}, points = {{-20, 100}, {0, 60}, {20, 100}}), Text(extent = {{-90, 100}, {-50, 80}}, textString = "%n_QForcing")}, coordinateSystem(initialScale = 0.1)));
+end QMForcing;

src/rtctools_channel_flow/modelica/Deltares/ChannelFlow/Salt/Elements/SaltyLinearReservoir.mo

@@ -0,0 +1,18 @@
+within Deltares.ChannelFlow.Salt.Elements;
+
+model SaltyLinearReservoir
+  /*
+  This block is designed to be used together with the "salt_simulation_mixin" to calculate dispersive and advective transport
+   between salty reservoir elements, do not user in optimization.
+  */
+  import SI = Modelica.Units.SI;
+  extends SaltyPartialReservoir(H(min = H_b));
+
+  parameter SI.Area A;
+  parameter SI.Position H_b; // Bed level
+
+equation
+
+  V = A * (H - H_b);
+  annotation(Icon(coordinateSystem(extent = {{-100, -100}, {100, 100}}, preserveAspectRatio = true, initialScale = 0.1, grid = {10, 10})));
+end SaltyLinearReservoir;

src/rtctools_channel_flow/modelica/Deltares/ChannelFlow/Salt/Elements/SaltyLinearReservoirBnd.mo

@@ -0,0 +1,18 @@
+within Deltares.ChannelFlow.Salt.Elements;
+
+model SaltyLinearReservoirBnd
+  /*
+  This block is designed to be used together with the "salt_simulation_mixin" to calculate dispersive and advective transport
+   between salty reservoir elements, do not user in optimization.
+  */
+  import SI = Modelica.Units.SI;
+  extends SaltyPartialReservoirBnd(H(min = H_b));
+
+  parameter SI.Area A;
+  parameter SI.Position H_b; // Bed level
+
+equation
+
+  V = A * (H - H_b);
+  annotation(Icon(coordinateSystem(extent = {{-100, -100}, {100, 100}}, preserveAspectRatio = true, initialScale = 0.1, grid = {10, 10})));
+end SaltyLinearReservoirBnd;

src/rtctools_channel_flow/modelica/Deltares/ChannelFlow/Salt/Elements/SaltyPartialReservoir.mo

@@ -0,0 +1,31 @@
+within Deltares.ChannelFlow.Salt.Elements;
+
+partial model SaltyPartialReservoir
+  /*
+  This block is designed to be used together with the "salt_simulation_mixin" to calculate dispersive and advective transport
+   between salty reservoir elements, do not user in optimization.
+  */
+  import SI = Modelica.Units.SI;
+  extends QMForcing;
+  extends Deltares.ChannelFlow.Internal.QLateral;
+  extends Deltares.ChannelFlow.Internal.Reservoir;
+
+  SI.Position H;
+  SI.Density C(min = 0, nominal = 1);
+
+  parameter SI.VolumeFlowRate Q_nominal = 1.0;
+  parameter SI.VolumeFlowRate C_nominal = 1.0;
+  parameter SI.Volume V_nominal;
+
+equation
+  H = HQUp.H;
+  H = HQDown.H;
+
+  C = HQUp.C[1];
+  C = HQDown.C[1];
+
+  // Mass balance, can be used for one medium - salt
+  der(V) = HQUp.Q + HQDown.Q + sum(QForcing) + sum(QLateral.Q);
+  der(V * C) = HQUp.M[1] + HQDown.M[1] + sum(MForcing[:,1]);
+
+end SaltyPartialReservoir;

src/rtctools_channel_flow/modelica/Deltares/ChannelFlow/Salt/Elements/SaltyPartialReservoirBnd.mo

@@ -0,0 +1,32 @@
+within Deltares.ChannelFlow.Salt.Elements;
+
+partial model SaltyPartialReservoirBnd
+  /*
+  This block is designed to be used together with the "salt_simulation_mixin" to calculate dispersive and advective transport
+   between salty reservoir elements, do not user in optimization.
+  */
+  import SI = Modelica.Units.SI;
+  extends Deltares.ChannelFlow.Internal.HQTwoPort;
+  extends Deltares.ChannelFlow.Internal.QForcing;
+  extends Deltares.ChannelFlow.Internal.QLateral;
+  extends Deltares.ChannelFlow.Internal.Reservoir;
+
+  SI.Position H;
+  SI.Density C[medium.n_substances](each min = 0, each nominal = 1);
+
+  parameter SI.VolumeFlowRate Q_nominal = 1.0;
+  parameter SI.VolumeFlowRate C_nominal = 1.0;
+  parameter SI.Volume V_nominal;
+
+equation
+
+  H = HQUp.H;
+  H = HQDown.H;
+
+  C = HQUp.C;
+  C = HQDown.C;
+
+  der(V) = 0;
+  der(V * C) = fill(0.0, medium.n_substances);
+
+end SaltyPartialReservoirBnd;

src/rtctools_channel_flow/modelica/Deltares/ChannelFlow/Salt/Elements/SaltyPartialStorage.mo

@@ -0,0 +1,22 @@
+within Deltares.ChannelFlow.Salt.Elements;
+
+partial model SaltyPartialStorage
+  /*
+  This block is designed to be used together with the "salt_simulation_mixin" to calculate dispersive and advective transport
+   between salty reservoir elements, do not user in optimization.
+  */
+  import SI = Modelica.Units.SI;
+  extends Deltares.ChannelFlow.Internal.HQOnePort(HQ.Q(each nominal = Q_nominal), HQ.M(each nominal =Q_nominal - C_nominal));
+  extends Deltares.ChannelFlow.Internal.QForcing(QForcing(each nominal = Q_nominal));
+  extends Deltares.ChannelFlow.Internal.Volume;
+
+  parameter SI.Volume V_nominal;
+  parameter SI.Density C_nominal = 1e-3;
+  parameter SI.VolumeFlowRate Q_nominal = 1.0;
+
+equation
+
+  der(V) / Q_nominal = (HQ.Q + sum(QForcing)) / Q_nominal;
+  HQ.M / (Q_nominal * C_nominal) =  der(V * HQ.C)   / (Q_nominal * C_nominal);
+
+end SaltyPartialStorage;

src/rtctools_channel_flow/modelica/Deltares/ChannelFlow/Salt/Elements/SubstanceControlledStructure.mo

@@ -0,0 +1,73 @@
+within Deltares.ChannelFlow.Salt.Elements;
+
+model SubstanceControlledStructure "SubstanceControlledStructure"
+  /*
+  This block is designed to be used together with the "salt_simulation_mixin" to calculate dispersive and advective transport
+   between salty reservoir elements, do not user in optimization.
+  */
+  extends Deltares.ChannelFlow.Internal.HQTwoPort;
+  import SI = Modelica.Units.SI;
+  function smooth_abs = Deltares.ChannelFlow.Internal.Functions.SmoothAbs;
+
+  parameter SI.MassFlowRate Q_nominal = 1;
+  parameter SI.Density C_nominal = 1e-3;
+  parameter SI.Height width = 2000;
+  parameter SI.Temperature temperature_up;
+  parameter SI.Temperature temperature_down;
+  parameter SI.Height H_b_up;
+  parameter SI.Height H_b_down;
+
+  SI.Density rho_up(nominal=1000, start = 1000.0);
+  SI.Density rho_down(nominal=1000, start = 1000.0);
+  SI.Concentration salinity_psu_up(nominal=34.7, start = 34.7);
+  SI.Concentration salinity_psu_down(nominal=34.7, start = 34.7);
+  SI.Density rho_ref_up;
+  SI.Density rho_ref_down;
+
+  Real a_up;
+  Real b_up;
+  Real c_up;
+  Real a_down;
+  Real b_down;
+  Real c_down;
+  Real flux_q1_s1;
+  Real epsilon_abs = 0.000001;
+
+equation
+
+  salinity_psu_up = HQUp.C[1] / rho_up * 1000.0;
+  salinity_psu_down = HQDown.C[1] / rho_down * 1000.0;
+
+  //Using UNESCO equation of state (EOS-80) provides a way to calculate the density of seawater as a function of salinity, temperature, and pressure.
+  a_up = 8.24493E-1 - 4.0899E-3 * temperature_up + 7.6438E-5 * temperature_up^2.0;// - 8.2467E-7 * temperature_up^3.0 + 5.3875E-9 * temperature_up^4.0;
+  b_up = -5.72466E-3 + 1.0227E-4 * temperature_up - 1.6546E-6 * temperature_up^2.0;
+  c_up = 4.8314E-4;
+  rho_ref_up = (999.842594 + 6.793952E-2 * temperature_up - 9.095290E-3 * temperature_up^2.0 +1.001685E-4 * temperature_up^3.0 - 1.120083E-6 *temperature_up^4.0 +6.536332E-9 * temperature_up^5.0);
+  rho_up = rho_ref_up + a_up * salinity_psu_up + b_up * salinity_psu_up^1.5 + c_up * salinity_psu_up^2.0;
+
+  a_down = 8.24493E-1 - 4.0899E-3 * temperature_down + 7.6438E-5 * temperature_down^2.0;// - 8.2467E-7 * temperature_down^3.0 + 5.3875E-9 * temperature_down^4.0;
+  b_down = -5.72466E-3 + 1.0227E-4 * temperature_down - 1.6546E-6 * temperature_down^2.0;
+  c_down = 4.8314E-4;
+  rho_ref_down = (999.842594 + 6.793952E-2 * temperature_down - 9.095290E-3 * temperature_down^2.0 +1.001685E-4 * temperature_down^3.0 - 1.120083E-6 *temperature_down^4.0 +6.536332E-9 * temperature_down^5.0); 
+  rho_down = rho_ref_down + a_down * salinity_psu_down + b_down * salinity_psu_down^1.5 + c_down * salinity_psu_down^2.0;
+
+  flux_q1_s1 =  (2*9.81)^0.5 * width / 2 * min(HQUp.H-H_b_up, HQDown.H-H_b_down)^1.5*(smooth_abs(rho_up-rho_down, epsilon_abs)/(rho_up+rho_down))^0.5;
+
+  /* This was a previous implementation, without allowing two-directional flow
+  if HQUp.Q  > flux_q1_s1 then
+      HQUp.M = HQUp.Q * HQUp.C[1];
+  else
+      HQUp.M =0.5 * HQUp.Q * (HQUp.C[1]+ HQDown.C[1]) + (HQUp.C[1]-HQDown.C[1])* 0.5 * (2*9.81)^0.5 * width / 2 * min(HQUp.H, HQDown.H)^1.5*(smooth_abs(rho_up-rho_down, epsilon_abs)/(rho_up+rho_down))^0.5;
+  end if;
+  */
+
+  if HQUp.Q < -flux_q1_s1 then
+     HQUp.M[1] = HQUp.Q * HQDown.C[1];
+  elseif HQUp.Q  > flux_q1_s1 then
+      HQUp.M[1] = HQUp.Q * HQUp.C[1];
+  else
+      HQUp.M[1] =0.5 * HQUp.Q * (HQUp.C[1]+ HQDown.C[1]) + (HQUp.C[1]-HQDown.C[1])* 0.5 * (2*9.81)^0.5 * width / 2 * min(HQUp.H-H_b_up, HQDown.H-H_b_down)^1.5*(smooth_abs(rho_up-rho_down, epsilon_abs)/(rho_up+rho_down))^0.5;
+  end if;
+
+  annotation(Icon(coordinateSystem( initialScale = 0.1, grid = {10, 10}), graphics = {Polygon(origin = {0, -16.67}, fillColor = {255, 128, 0}, fillPattern = FillPattern.Solid, points = {{0, 66.667}, {-50, -33.333}, {50, -33.333}, {0, 66.667}})}), Diagram(coordinateSystem(extent = {{-100, -100}, {100, 100}}, preserveAspectRatio = true, initialScale = 0.1, grid = {10, 10})));
+end SubstanceControlledStructure;

src/rtctools_channel_flow/modelica/Deltares/ChannelFlow/Salt/Elements/package.mo

@@ -0,0 +1,4 @@
+within Deltares.ChannelFlow.Salt;
+
+package Elements
+end Elements;

src/rtctools_channel_flow/modelica/Deltares/ChannelFlow/Salt/Elements/package.order

@@ -0,0 +1,8 @@
+SubstanceControlledStructure
+SaltyPartialStorage
+SaltyPartialReservoirBnd
+SaltyPartialReservoir
+SaltyLinearReservoirBnd
+SaltyLinearReservoir
+QMForcing
+NodeSalty

src/rtctools_channel_flow/modelica/Deltares/ChannelFlow/Salt/package.mo

@@ -0,0 +1,4 @@
+within Deltares.ChannelFlow;
+
+package Salt
+end Salt;

src/rtctools_channel_flow/modelica/Deltares/ChannelFlow/Salt/package.order

@@ -0,0 +1 @@
+Elements

src/rtctools_channel_flow/modelica/Deltares/ChannelFlow/package.order

@@ -3,3 +3,4 @@ Interfaces
 Internal
 SimpleRouting
 Hydraulic
+Salt