Started adding new door model

For #1353 [ci skip]
ibpsa · Oct 7, 2020 · f54dd33 · f54dd33
1 parent b17ef04
commit f54dd33
Show file tree

Hide file tree

Showing 2 changed files with 313 additions and 0 deletions.
diff --git a/IBPSA/Airflow/Multizone/BaseClasses/Door.mo b/IBPSA/Airflow/Multizone/BaseClasses/Door.mo
@@ -0,0 +1,128 @@
+within IBPSA.Airflow.Multizone.BaseClasses;
+partial model Door
+  "Door model for bi-directional flow"
+  extends IBPSA.Fluid.Interfaces.PartialFourPortInterface(
+    redeclare final package Medium1 = Medium,
+    redeclare final package Medium2 = Medium,
+    final allowFlowReversal1=true,
+    final allowFlowReversal2=true,
+    final m1_flow_nominal=10/3600*1.2,
+    final m2_flow_nominal=m1_flow_nominal,
+    final m1_flow_small=1E-4*abs(m1_flow_nominal),
+    final m2_flow_small=1E-4*abs(m2_flow_nominal));
+
+  replaceable package Medium =
+    Modelica.Media.Interfaces.PartialMedium "Medium in the component"
+      annotation (choices(
+        choice(redeclare package Medium = IBPSA.Media.Air "Moist air")));
+
+  parameter Modelica.SIunits.Length wOpe=0.9 "Width of opening"
+    annotation (Dialog(group="Geometry"));
+  parameter Modelica.SIunits.Length hOpe=2.1 "Height of opening"
+    annotation (Dialog(group="Geometry"));
+
+  parameter Modelica.SIunits.PressureDifference dp_turbulent(
+    min=0,
+    displayUnit="Pa") = 0.01
+    "Pressure difference where laminar and turbulent flow relation coincide. Recommended: 0.01";
+
+  Modelica.SIunits.VolumeFlowRate VAB_flow(nominal=0.001)
+    "Volume flow rate from A to B if positive";
+  Modelica.SIunits.VolumeFlowRate VBA_flow(nominal=0.001)
+    "Volume flow rate from B to A if positive";
+  Modelica.SIunits.MassFlowRate mAB_flow(nominal=0.001)
+    "Mass flow rate from A to B if positive";
+  Modelica.SIunits.MassFlowRate mBA_flow(nominal=0.001)
+    "Mass flow rate from B to A if positive";
+
+  Modelica.SIunits.Velocity vAB(nominal=0.01) "Average velocity from A to B";
+  Modelica.SIunits.Velocity vBA(nominal=0.01) "Average velocity from B to A";
+
+protected
+  parameter Medium.ThermodynamicState sta_default=Medium.setState_pTX(
+      T=Medium.T_default,
+      p=Medium.p_default,
+      X=Medium.X_default);
+
+  parameter Modelica.SIunits.Density rho_default=Medium.density(sta_default)
+    "Density";
+
+  Modelica.SIunits.Density rho_a1_inflow
+    "Density of air flowing in from port_a1";
+  Modelica.SIunits.Density rho_a2_inflow
+    "Density of air flowing in from port_a2";
+
+equation
+  // Compute the density of the inflowing media.
+  rho_a1_inflow = Medium1.density(state_a1_inflow);
+  rho_a2_inflow = Medium2.density(state_a2_inflow);
+
+  mAB_flow = rho_a1_inflow*VAB_flow;
+  mBA_flow = rho_a2_inflow*VBA_flow;
+  // Average velocity (using the whole orifice area)
+  vAB = VAB_flow/A;
+  vBA = VBA_flow/A;
+
+  port_a1.m_flow = mAB_flow;
+  port_a2.m_flow = mBA_flow;
+
+  // Energy balance (no storage, no heat loss/gain)
+  port_a1.h_outflow = inStream(port_b1.h_outflow);
+  port_b1.h_outflow = inStream(port_a1.h_outflow);
+  port_a2.h_outflow = inStream(port_b2.h_outflow);
+  port_b2.h_outflow = inStream(port_a2.h_outflow);
+
+  // Mass balance (no storage)
+  port_a1.m_flow = -port_b1.m_flow;
+  port_a2.m_flow = -port_b2.m_flow;
+
+  port_a1.Xi_outflow = inStream(port_b1.Xi_outflow);
+  port_b1.Xi_outflow = inStream(port_a1.Xi_outflow);
+  port_a2.Xi_outflow = inStream(port_b2.Xi_outflow);
+  port_b2.Xi_outflow = inStream(port_a2.Xi_outflow);
+
+  // Transport of trace substances
+  port_a1.C_outflow = inStream(port_b1.C_outflow);
+  port_b1.C_outflow = inStream(port_a1.C_outflow);
+  port_a2.C_outflow = inStream(port_b2.C_outflow);
+  port_b2.C_outflow = inStream(port_a2.C_outflow);
+
+  annotation (
+    Icon(graphics={
+        Rectangle(
+          extent={{-60,80},{60,-84}},
+          lineColor={0,0,255},
+          pattern=LinePattern.None,
+          fillColor={85,75,55},
+          fillPattern=FillPattern.Solid),
+        Rectangle(
+          extent={{-54,72},{56,-84}},
+          lineColor={0,0,0},
+          fillColor={215,215,215},
+          fillPattern=FillPattern.Solid),
+        Polygon(
+          points={{56,72},{-36,66},{-36,-90},{56,-84},{56,72}},
+          lineColor={0,0,0},
+          fillColor={95,95,95},
+          fillPattern=FillPattern.Solid),
+        Polygon(
+          points={{-30,-10},{-16,-8},{-16,-14},{-30,-16},{-30,-10}},
+          lineColor={0,0,255},
+          pattern=LinePattern.None,
+          fillColor={0,0,0},
+          fillPattern=FillPattern.Solid)}),
+    Documentation(info="<html>
+<p>
+This is a partial model for the bi-directional air flow through a door.
+</p>
+</html>",
+revisions="<html>
+<ul>
+<li>
+October 6, 2020, by Michael Wetter:<br/>
+First implementation for
+<a href=\"https://github.com/ibpsa/modelica-ibpsa/issues/1353\">#1353</a>.
+</li>
+</ul>
+</html>"));
+end Door;
diff --git a/IBPSA/Airflow/Multizone/DoorOpen.mo b/IBPSA/Airflow/Multizone/DoorOpen.mo
@@ -0,0 +1,185 @@
+within IBPSA.Airflow.Multizone;
+model DoorOpen
+  "Door model using discretization along height coordinate"
+  extends IBPSA.Airflow.Multizone.BaseClasses.Door;
+
+  parameter Real CD=0.65 "Discharge coefficient"
+    annotation (Dialog(group="Orifice characteristics"));
+
+  final Modelica.SIunits.Area A = wOpe*hOpe "Face area";
+
+protected
+  constant Real mFixed = 0.5 "Fixed value for flow coefficient";
+  constant Real gamma(min=1) = 1.5
+    "Normalized flow rate where dphi(0)/dpi intersects phi(1)";
+  constant Real a = gamma
+    "Polynomial coefficient for regularized implementation of flow resistance";
+  constant Real b = 1/8*mFixed^2 - 3*gamma - 3/2*mFixed + 35.0/8
+    "Polynomial coefficient for regularized implementation of flow resistance";
+  constant Real c = -1/4*mFixed^2 + 3*gamma + 5/2*mFixed - 21.0/4
+    "Polynomial coefficient for regularized implementation of flow resistance";
+  constant Real d = 1/8*mFixed^2 - gamma - mFixed + 15.0/8
+    "Polynomial coefficient for regularized implementation of flow resistance";
+
+  parameter Real kVal=CD*A*sqrt(2/rho_default) "Flow coefficient, k = V_flow/ dp^m";
+  parameter Reak kT = CD*wOpe*sqrt(2/rho_default)
+    *(g*rho_default*hOpe/2/Medium.T_default)^mFixed *hOpe/2
+    / conTP^mFixed
+    "Constant coefficient for buoyancy driven air flow rate";
+  parameter Real conTP = IBPSA.Media.Air.dStp*Modelica.Media.IdealGases.Common.SingleGasesData.Air.R
+    "Conversion factor for converting temperature difference to pressure difference";
+
+
+equation
+  // Air flow rate due to static pressure difference
+  VABp_flow = IBPSA.Airflow.Multizone.BaseClasses.powerLawFixedM(
+      k=kVal,
+      dp=port_a1.p-port_a2.p,
+      m=mFixed,
+      a=a,
+      b=b,
+      c=c,
+      d=d,
+      dp_turbulent=dp_turbulent);
+  // Air flow rate due to buoyancy
+  // Because powerLawFixedM requires as an input a pressure difference pa-pb,
+  // we convert Ta-Tb by multiplying it with rho*R, and we divide
+  // above the constant expression by (rho*R)^m on the right hand-side of kT.
+  VABt_flow = IBPSA.Airflow.Multizone.BaseClasses.powerLawFixedM(
+      k=kT,
+      dp=conTP*(Medium.temperature(state_a1_inflow)-Medium.temperature(state_a2_inflow)),
+      m=mFixed,
+      a=a,
+      b=b,
+      c=c,
+      d=d,
+      dp_turbulent=dp_turbulent);
+  // Net flow rate
+  // fixme: check sign and make sure documentation is consistent.
+  VAB_flow = +VABp_flow/2 + VABt_flow;
+  VBA_flow = -VABp_flow/2 - VABt_flow;
+
+  annotation (defaultComponentName="doo",
+Documentation(info="<html>
+<p>
+Model for bi-directional air flow through a large opening such as a door.
+</p>
+<p>
+The air flow is composed of two components, a one-directional bulk air flow
+due to static pressure difference in the adjoining two thermal zones, and
+a two-directional airflow due to temperature-induced differences in density
+of the air in the two rooms.
+Although turbulent air flow is a nonlinear phenomenon and hence the
+two air flow rates could only be superposed for laminar flow,
+the model is based on the simplifying assumption that these two
+air flow rates can be superposed.
+This assumption is made because
+it leads to a simple model and because there is significant uncertainty
+and assumptions anyway in the model for bidirectional flow through a door.
+</p>
+<h4>Main equations</h4>
+<p>
+The air flow rate due to static pressure difference is
+</p>
+<p align=\"center\" style=\"font-style:italic;\">
+    V&#775;<sub>ab,p</sub> =  C<sub>D</sub> w h (2/&rho;<sub>0</sub>)<sup>0.5</sup>  &Delta;p<sup>m</sup>,
+</p>
+<p>
+where
+<i>V&#775;</i> is the volume flow rate,
+<i>C<sub>D</sub></i> is the discharge coefficient,
+<i>w</i> and <i>h</i> are the width and height of the opening
+<i>&rho;<sub>0</sub></i> is the mass density at the medium default pressure, temperature and humidity,
+<i>m</i> is the flow exponent and
+<i>&Delta;p = p<sub>a</sub> - p<sub>b</sub></i> is the static pressure difference between
+the thermal zones.
+For this discussion, we will assume <i>p<sub>a</sub> &gt; p<sub>b</sub></i>.
+For turbulent flow, <i>m=1/2</i> and for laminar flow <i>m=1</i>.
+</p>
+<p>
+The air flow rate due to temperature difference in the rooms is
+<i>V&#775;<sub>ab,t</sub></i> for flow from thermal zone <i>a</i> to <i>b</i>,
+and
+<i>V&#775;<sub>ba,t</sub></i> for flow from thermal zone <i>b</i> to <i>a</i>.
+The net air flow rate from <i>a</i> to <i>b</i> is
+</p>
+<p align=\"center\" style=\"font-style:italic;\">
+    V&#775;<sub>ab</sub> = +&nbsp; V&#775;<sub>ab,p</sub>/2 + &nbsp; V&#775;<sub>ab,t</sub>,
+</p>
+and, similarly,
+<p align=\"center\" style=\"font-style:italic;\">
+    V&#775;<sub>ba</sub> = -&nbsp; V&#775;<sub>ab,p</sub>/2 + &nbsp; V&#775;<sub>ba,t</sub>.
+</p>
+<p>
+To calculate <i>V&#775;<sub>ba,p</sub></i>, we assume
+</p>
+<p align=\"center\" style=\"font-style:italic;\">
+    m&#775;<sub>ab,t</sub> =  -m&#775;<sub>ba,t</sub> = &rho;<sub>0</sub> &nbsp; V&#775;<sub>ab,t</sub> =  -&rho;<sub>0</sub> &nbsp; V&#775;<sub>ba,t</sub>.
+</p>
+<p>
+With this assumption, the neutral height, e.g., the height where the air flow rate due to flow
+induced by temperature difference is zero, is at <i>h/2</i>.
+Hence,
+</p>
+<p align=\"center\" style=\"font-style:italic;\">
+    V&#775;<sub>ab,t</sub> =  &int;<sub>0</sub><sup>h/2</sup> (&part; &frasl; &part; z) V&#775;<sub>ab,t</sub>(z) &nbsp; dz,
+</p>
+<p>
+and with
+</p>
+<p align=\"center\" style=\"font-style:italic;\">
+    V&#775;<sub>ab,t</sub> =  C<sub>D</sub> w h (2/&rho;<sub>0</sub>)<sup>0.5</sup>  &Delta;p<sub>ab</sub><sup>m</sup>(z),
+</p>
+and
+<p align=\"center\" style=\"font-style:italic;\">
+  &Delta;p<sub>ab</sub><sup>m</sup>(z) = g z (&rho;<sub>a</sub>-&rho;<sub>b</sub>) = g z (&rho;<sub>a</sub>-&rho;<sub>b</sub>)
+    &approx; &rho;<sub>0</sub> (T<sub>b</sub> - T<sub>a</sub>) &frasl; T<sub>0</sub>,
+</p>
+<p>
+where we used
+<i>&rho;<sub>a</sub> = p<sub>0</sub> /(R T<sub>a</sub>)</i> and
+<i>T<sub>a</sub> T<sub>b</sub> &approx; T<sub>0</sub></i>.
+Substituting this expression into the integral yields
+</p>
+<p align=\"center\" style=\"font-style:italic;\">
+    V&#775;<sub>ab,t</sub> =  C<sub>D</sub> w (2/&rho;<sub>0</sub>)<sup>0.5</sup> g<sup>m</sup> &rho;<sub>0</sub><sup>m</sup>
+      (h/2)<sup>m+1</sup> (T<sub>b</sub>-T<sub>a</sub>)<sup>m</sup> &frasl; T<sub>0</sub><sup>m</sup>.
+</p>
+<h4>Main assumptions</h4>
+<p>
+The main assumptions are
+</p>
+<ul>
+<li>
+<p>
+that the air flow rates due to static pressure difference and due to temperature-difference can be superposed,
+</p>
+</li>
+<li>
+<p>
+that for computing the neutral height, <i>m&#775;<sub>ab,t</sub> =  &rho;<sub>0</sub> &nbsp; V&#775;<sub>ab,t</sub></i>, and
+</p>
+</li>
+<li>
+<p>
+that for buoyancy-driven flow the air density can be approximated as <i>&rho; = p<sub>0</sub> /(R T)</i>,
+</p>
+</li>
+</ul>
+<h4>Notes</h4>
+<p>
+For a more detailed model, use
+<a href=\"modelica://IBPSA.Airflow.Multizone.DoorDiscretizedOpen\">
+IBPSA.Airflow.Multizone.DoorDiscretizedOpen</a>.
+</p>
+</html>",
+revisions="<html>
+<ul>
+<li>
+October 6, 2020, by Michael Wetter:<br/>
+First implementation for
+<a href=\"https://github.com/ibpsa/modelica-ibpsa/issues/1353\">#1353</a>.
+</li>
+</ul>
+</html>"));
+end DoorOpen;