diff --git a/FluidDynamics/Components/Orifices/OpeningAnalytic.mo b/FluidDynamics/Components/Orifices/OpeningAnalytic.mo
new file mode 100644
index 0000000000000000000000000000000000000000..9ac507891b62047261c79b7f28f033e43314edeb
--- /dev/null
+++ b/FluidDynamics/Components/Orifices/OpeningAnalytic.mo
@@ -0,0 +1,99 @@
+within TAeZoSysPro.FluidDynamics.Components.Orifices;
+
+model OpeningAnalytic
+  replaceable package Medium = TAeZoSysPro.Media.MyMedia ;
+  
+  // User defined parameters
+  parameter Real Cd = 0.61 "discharge coefficient";
+  parameter Modelica.SIunits.CrossSection A = 1 "Opening cross section";
+  parameter Modelica.SIunits.Height H = 1 "Opening's Height";   
+  parameter Modelica.SIunits.Length Alt_a = 1.0 "Altitude of port_a";
+  parameter Modelica.SIunits.Length Alt_b = 1.0 "Altitude of port_b";
+  parameter Modelica.SIunits.Length Alt_opening = 1.0 "Altitude at the opening center";
+    
+  // Internal variables
+  Modelica.SIunits.Pressure p_a "Pressure of node connected to port_a";
+  Modelica.SIunits.Pressure p_b "Pressure of node connected to port_b";
+  Modelica.SIunits.PressureDifference dp;   
+  Modelica.SIunits.Density rho_A,rho_B;
+  Medium.Density[Medium.nX] rho_up, rho_down "Upstream density in upper and lower part of the opening";
+  Medium.MassFraction[Medium.nX] X_up, X_down;
+  Modelica.SIunits.Height HN "Heigh of the point where flow switches in direction (zero flow)";
+  Medium.ThermodynamicState state_a, state_b;
+  Modelica.SIunits.MassFlowRate m_flow_up, m_flow_down "mass flow rate in upper and lower part of the opening";
+  
+  // Imported modules
+  TAeZoSysPro.FluidDynamics.Interfaces.FlowPort_a port_a(redeclare package Medium = Medium) annotation (
+    Placement(visible = true, transformation(origin = {-58, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {-90, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
+  TAeZoSysPro.FluidDynamics.Interfaces.FlowPort_b port_b(redeclare package Medium = Medium) annotation (
+    Placement(visible = true, transformation(origin = {38, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {90, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
+
+protected
+  Modelica.SIunits.SpecificEnthalpy h_a "Specific enthalpy from port_a" ;
+  Modelica.SIunits.SpecificEnthalpy h_b "Specific enthalpy from port_b" ;
+  parameter Modelica.SIunits.Velocity Vel_small = 0.001 ;
+  constant Modelica.SIunits.Acceleration g = Modelica.Constants.g_n "Gravitational acceleration";  
+  
+equation
+//
+  state_a = Medium.setState_dTX(d = sum(port_a.d), T = port_a.T, X = port_a.d / sum(port_a.d));
+  state_b = Medium.setState_dTX(d = sum(port_b.d), T = port_b.T, X = port_b.d / sum(port_b.d));
+
+// pressure reconstruction
+  p_a = Medium.pressure(state_a) + sum(port_a.d) * g * (Alt_a-Alt_opening);
+  p_b = Medium.pressure(state_b) + sum(port_b.d) * g * (Alt_b-Alt_opening);
+  dp = p_a - p_b;
+  
+// specific enthalpy reconstruction
+  h_a = Medium.specificEnthalpy(state_a);
+  h_b = Medium.specificEnthalpy(state_b);
+  
+// density
+  rho_A = sum(port_a.d);
+  rho_B = sum(port_b.d);
+
+// upsteam density
+  rho_up = smooth(1, if dp - (sum(port_a.d) - sum(port_b.d)) * g * (H - HN) > 0 then port_a.d else port_b.d);
+  rho_down = port_a.d + port_b.d - rho_up;
+  
+// upstream mass fraction
+  X_up = rho_up / sum(rho_up);
+  X_down = port_a.d/rho_A + port_b.d/rho_B - X_up;
+  
+// altidude of change in direction of flow (if change happens)
+  HN = min(max(dp / ((rho_A - rho_B) * g) + H/2, 0), H);
+
+// m_flow = f(dp, Δρ, ΔH)
+  m_flow_up = Cd * A/H * sqrt(2 * sum(rho_up)) * ( Modelica.Fluid.Utilities.regRoot(x=dp, delta=1e-3) * (H - HN) - Modelica.Fluid.Utilities.regRoot(x=(rho_A - rho_B) * g, delta=1e-3) * 2/3 * (H-HN)^(3/2) );
+  m_flow_down = Cd * A/H * sqrt(2 * sum(rho_down)) * ( Modelica.Fluid.Utilities.regRoot(x=dp, delta=1e-3) * (- HN) - Modelica.Fluid.Utilities.regRoot(x=(rho_A - rho_B) * g, delta=1e-3) * 2/3 * (-HN^(3/2)) );
+  
+// Port handover
+  port_a.m_flow = m_flow_up * X_up + m_flow_down * X_down ;
+  port_a.m_flow + port_b.m_flow = fill(0.0, Medium.nX);
+  
+  port_a.H_flow = semiLinear(m_flow_up, h_a + g * (Alt_a - Alt_b), h_b - g * (Alt_a - Alt_b));
+  port_a.H_flow + port_b.H_flow = 0;
+  
+  annotation(
+    Icon(graphics = {Line(origin = {0, 69.8886}, points = {{0, 30}, {0, -30}}, thickness = 2), Line(origin = {0, -70}, points = {{0, 30}, {0, -30}}, thickness = 2)}, coordinateSystem(initialScale = 0.1)),
+    experiment(StartTime = 0, StopTime = 1000, Tolerance = 1e-06, Interval = 1),
+    Documentation(info = "<html>
+  <head>
+    <title>Opening</title>
+  </head>
+	
+  <body lang=\"en-UK\">
+    <p>
+      This components allows to model the mass flow rate through from either static boundary pressure difference or buoyancy effect through a vertical orifice in a wall spliting two ambiances.
+    </p>
+    
+    <p>
+      To be considered as an orifice, the depth of the hole in the wall has to remain bellow the hydrodynamic entrance region (Distance between the entrance of the hole and the position where the dynamic boundary layers meet).
+    </p>
+
+    <p>    
+      <a href=\"modelica://TAeZoSysPro/Information/FluidDynamics/Components/Orifices/build/Demo_opening_aerodynamic.pdf\">See demonstration</a>.
+    </p>   
+</html>"));
+
+end OpeningAnalytic;
\ No newline at end of file
diff --git a/FluidDynamics/Components/Orifices/package.order b/FluidDynamics/Components/Orifices/package.order
index 8f428a126ab065ef4d8d3e0e3b048979ffd9f83b..6703c422b9b993b5042dc4ec9e36973dc03294af 100644
--- a/FluidDynamics/Components/Orifices/package.order
+++ b/FluidDynamics/Components/Orifices/package.order
@@ -1,3 +1,4 @@
 Opening
 HorizontalOpening
 SimpleOpeningComp
+OpeningAnalytic
diff --git a/Information/FluidDynamics/Components/Orifices/Demo_opening_aerodynamic.tex b/Information/FluidDynamics/Components/Orifices/Demo_opening_aerodynamic.tex
new file mode 100644
index 0000000000000000000000000000000000000000..dcdbb3c948659226aba9004f055336f626742329
--- /dev/null
+++ b/Information/FluidDynamics/Components/Orifices/Demo_opening_aerodynamic.tex
@@ -0,0 +1,132 @@
+\documentclass[a4paper, french]{article}
+\usepackage[utf8]{inputenc}
+\usepackage[T1]{fontenc}
+
+% package pour gérer la langue du document
+\usepackage{babel}
+
+
+\usepackage{amsmath,amsfonts,amssymb}
+
+% Titre du document
+\title{Vertical Opening Thermo Aeraulic}
+
+\begin{document}
+    \maketitle         
+            
+    \section{$0 \leq z \leq H $}
+    
+        \subsection{Mass flow rate}
+            
+            To quantity the mass flow rates of transfer, first express the diffential of pressure over the height.
+
+            \begin{equation}
+                \label{eq1}
+                \mathrm{d}p = -\rho \times g \times \mathrm{d} z
+            \end{equation}
+            
+            Now establish the pressure difference of the gas colomn from the two ambiances A and B.
+            
+            \begin{equation}
+                \label{eq2}
+                \mathrm{d}p_{AB} = \Delta {p_{AB}}_{static} -(\rho_A - \rho_B) \times g \times \mathrm{d} z
+            \end{equation}       
+            
+            
+            The altitude of the opening range from $ 0 \leq z \leq H $. The altitude where flow changes in direction is defined as $HN$. The integration of \eqref{eq2} is performed from $ z = HN $ to the bottom at $ z = 0 $
+
+            
+            \begin{equation}
+                \label{eq3}
+                \Delta p_{AB}(z) = \Delta {p_{AB}}_{static} + \int_{HN}^{z} -(\rho_A - \rho_B) \times g \times dx = -(\rho_A - \rho_B) \times g \times (z - HN)
+            \end{equation} 
+            
+            Considering that all the pressure difference is converted to dynamic pressure.
+            
+            \begin{equation}
+                \label{eq4}
+                \Delta p_{AB}(z) = \frac{1}{2} \times \rho \times Vel(z)^2
+            \end{equation}
+            
+            Now a problem arises with the choice of the density to use. The upstream density is used but can differs between the top and bottom part of the opening if the static pressure is lower than the hydrostatique pressure. Thefore the opening is cut in two and the altitude to which the flow changes in direction is called $ HN $.
+            
+            The relation between the flow velocity and the altitude derives:
+            
+            \begin{equation}
+                \label{eq5}
+                Vel(z)^2 = \left | \frac{2 \times \Delta {p_{AB}}_{static} - 2 \times (\rho_A - \rho_B) \times g \times (z-HN)}{\rho_{upstream} } \right |
+            \end{equation}
+            
+            The elementary mass flow rate derives from the velocity and the upstream density. The width of the opening is defined as $ w = \frac{A}{H} $:
+            
+            \begin{equation}
+                \label{eq6}
+                \mathrm{d} m\_flow = \rho_{upstream} \times Vel(z) \times Cd \times \frac{A}{H} \times \mathrm{d} z
+            \end{equation}
+            
+            Combining \eqref{eq5} and \eqref{eq6} to obtain:
+            
+            
+            \begin{equation}
+                \label{eq7}
+                \mathrm{d} m\_flow = Cd \times \frac{A}{H} \times \sqrt{ 2 \times \rho_{upstream} } \times \sqrt{\Delta {p_{AB}}_{static} - ( \rho_A - \rho_B ) \times g \times (z - HN)} \times \mathrm{d} z
+            \end{equation}
+            
+            Intregration of \eqref{eq7} over the height
+            
+            \begin{equation}
+                \label{eq8}
+                m\_flow = \int_{HN}^{z} \mathrm{d} m\_flow
+            \end{equation}
+            
+            From the height $ HN $ to  $ 0 $  (bottom part of the opening)
+            
+            \begin{equation} 
+            \label{eq9}
+                m\_flow_{down} = Cd \times \frac{A}{H} \times \sqrt{ 2 \times \rho_{down}} \times \left( \sqrt{ \Delta {p_{AB}} } \times (- HN) - \sqrt{ ( \rho_A - \rho_B ) \times g } \times \frac{2}{3} \times {(-HN)}^\frac{3}{2} \right)
+            \end{equation} 
+            
+            From the height $ HN $ to $ H $  (top part of the opening)
+            
+            \begin{equation} 
+            \label{eq10}
+                m\_flow_{up} = Cd \times \frac{A}{H} \times \sqrt{ 2 \times \rho_{up}} \times \left( \sqrt{ \Delta {p_{AB}} } \times (H-HN) - \sqrt{ ( \rho_A - \rho_B ) \times g } \times \frac{2}{3} \times {( H - HN)}^\frac{3}{2} \right)
+            \end{equation}     
+        
+            A modified version of the square root is used to deal with a negative value for the argument.
+            
+
+        \subsection{Zero pressure difference altitude calculation}
+    
+            Demonstration of the calculation of the altitude $ HN $ that leads to $ \Delta p_{AB}(HN) = 0 $ for a vertical opening. $ HN $ is limited to $ 0 \leq HN \leq H $.
+            The formula bellow reminds the mass flow rate at the upper part of the opening.
+            
+            \begin{equation*} 
+                m\_flow_{up} = Cd \times \frac{A}{H} \times \sqrt{ 2 \times \rho_{up}} \times \left( \sqrt{ \Delta {p_{AB}} } \times (H-HN) - \sqrt{ ( \rho_A - \rho_B ) \times g } \times \frac{2}{3} \times {( H - HN)}^\frac{3}{2} \right)
+            \end{equation*} 
+
+            The formula bellow reminds the mass flow rate at the lower part of the opening.
+
+            \begin{equation*} 
+                m\_flow_{down} = Cd \times \frac{A}{H} \times \sqrt{ 2 \times \rho_{down}} \times \left( \sqrt{ \Delta {p_{AB}} } \times (- HN) - \sqrt{ ( \rho_A - \rho_B ) \times g } \times \frac{2}{3} \times {(-HN)}^\frac{3}{2} \right)
+            \end{equation*} 
+            
+            Calculation of the height $ HN $ for the pressure difference between the opening to be zero from \eqref{eq3}. If there is no pressure difference, the altitude $ HN $ is $ \frac{H}{2} $. Therefore the lower bound of the integral is $ \frac{H}{2} $.
+            
+            \begin{equation*}
+                \Delta {p_{AB}}_{static} = \int_{\frac{H}{2}}^{HN} -(\rho_A - \rho_B) \times g \times dx = -(\rho_A - \rho_B) \times g \times (HN - \frac{H}{2})
+            \end{equation*}
+
+            \begin{equation}
+                HN = \frac{ \Delta {p_{AB}}_{static} }{ (\rho_A - \rho_B) \times g } + \frac{H}{2}
+            \end{equation}
+            
+            $HN$ is bounded to $ 0 \leq HN \leq H$, therefore:
+            
+            \begin{equation}
+                HN = min \left( max \left( \frac{ \Delta {p_{AB}}_{static} }{ (\rho_A - \rho_B) \times g } + \frac{H}{2}, 0\right), H \right)
+            \end{equation}            
+            
+        
+    
+\end{document}