Merge pull request #37 from andr1976/Valve-equations

Prep for in-house UV calc to allow HEOS mixtures

examples/ng.yml

@@ -0,0 +1,26 @@
+vessel:
+  length: 5.5
+  diameter: 0.344
+  thickness: 0.031
+  heat_capacity: 500
+  density: 7800. 
+  orientation: "horizontal"
+initial:
+  temperature: 273.0
+  pressure: 15500000.
+  fluid: "CH4" 
+calculation:
+  type: "energybalance"
+  time_step: 0.05
+  end_time: 300.
+valve:
+  flow: "discharge"
+  type: "orifice"
+  diameter: 0.009739
+  discharge_coef: 0.84
+  back_pressure: 101300.
+heat_transfer:
+  type: "specified_h"
+  temp_ambient: 273.
+  h_outer: 5 
+  h_inner: 'calc' 

examples/psv_sb.yml

@@ -4,17 +4,15 @@ vessel:
   orientation: "vertical"
   heat_capacity: 500
   density: 7700
-  orientation: "horizontal"
   thickness: 0.010 
 initial:
   temperature: 298.15
   pressure: 1601000.
   fluid: "N2" 
 calculation:
   type: "energybalance"
-  time_step: 5.
+  time_step: 1
   end_time: 1500.
-  eta: 1
 valve:
   flow: "discharge"
   type: "psv"

hyddown/hdclass.py

@@ -4,6 +4,8 @@
 
 import math
 import numpy as np
+from scipy.optimize import fmin
+from scipy.optimize import minimize
 from CoolProp.CoolProp import PropsSI
 from hyddown import transport as tp
 from hyddown import validator 
@@ -21,7 +23,7 @@ def __init__(self, input):
 
     def validate_input(self):
         valid = validator.validation(self.input)
-        if valid == False:
+        if valid is False:
             raise ValueError("Error in input file")
 
 
@@ -140,6 +142,25 @@ def initialize(self):
         self.m0 = self.rho0 * self.vol
 
 
+    def UDproblem(self, U, rho, Pguess, Tguess):
+        if "&" in self.species:                     
+            res_fun = lambda x: ((U-PropsSI("U","P",x[0],"T",x[1],self.species))/U)**2 + ((rho-PropsSI("D","P",x[0],"T",x[1],self.species))/rho)**2
+            x0=[Pguess,Tguess]
+            bounds=[(Pguess*0.95,Pguess+1e5),(Tguess-1,Tguess+1)]
+            res=minimize(res_fun,x0,method='Nelder-Mead')
+            P1 = res.x[0]
+            T1 = res.x[1]
+        else:
+            P1 = PropsSI(
+                "P", "D", rho, "U", U, self.species
+                )
+            T1 = PropsSI(
+                "T", "D", rho, "U", U, self.species
+               )
+
+        return P1, T1
+
+
     def run(self):
         # Inititialise
         input = self.input
@@ -274,11 +295,6 @@ def run(self):
                     "P", "D", self.rho[i], "U", self.U_mass[i - 1], self.species
                 )
             elif self.method == "energybalance":
-                P1 = PropsSI(
-                    "P", "D", self.rho[i], "T", self.T_fluid[i - 1], self.species
-                )
-                T1 = PropsSI("T", "P", P1, "H", self.H_mass[i - 1], self.species)
-
                 if self.heat_method == "specified_h" or self.heat_method == "detailed":
                     if self.h_in == "calc":
                         if self.vessel_orientation == "horizontal":
@@ -376,11 +392,6 @@ def run(self):
                 h_in = x * PropsSI("H", "T", T_used, "P", P_used, self.species) + (
                     1 - x
                 ) * PropsSI("U", "T", T_used, "P", P_used, self.species)
-                T_in = PropsSI("T", "P", self.P[i - 1], "H", h_in, self.species)
-                D_in = PropsSI("D", "P", self.P[i - 1], "T", T_in, self.species)
-                v_in = (self.mass_rate[i - 1] / D_in) / (
-                    self.diameter ** 2 / 4 * math.pi
-                )
                 P2 = self.P[i - 1]
                 if i > 1:
                     P1 = self.P[i - 2]
@@ -390,45 +401,11 @@ def run(self):
                     U_start
                     - self.tstep * self.mass_rate[i - 1] * h_in
                     + self.tstep * self.Q_inner[i]
-                )  # + (P2-P1) * self.vol
-                # print(U_start/self.mass_fluid[i-1],U_end/self.mass_fluid[i],h_in*self.mass_rate[i],self.Q_inner[i])
+                )  
                 self.U_mass[i] = U_end / self.mass_fluid[i]
-                P1 = PropsSI(
-                    "P", "D", self.rho[i], "U", U_end / self.mass_fluid[i], self.species
-                )
-                T1 = PropsSI(
-                    "T", "D", self.rho[i], "U", U_end / self.mass_fluid[i], self.species
-                )
-                # print('P1',P1,'T1',T1)
-                # U_start = NMOL_ADD * PropsSI('HMOLAR', 'P', P_used, 'T', T_used, self.species) + NMOL * PropsSI('HMOLAR', 'P', self.P[i-1], 'T', self.T_fluid[i-1], self.species) - self.eta *  self.P[i-1] * self.vol + self.Q_inner[i] * self.tstep
-                # NMOL=NMOL+NMOL_ADD
-
-                # U = 0
-                # nn = 0
-                # rho1 = 0
-                # itermax = 1000
-                # m = 0
-                # n = 0
-                # relax = 0.1
-
-                # while abs(self.rho[i] - rho1) > 0.01 and m < itermax:
-                #     m = m + 1
-                #     rho1 = PropsSI('D', 'T', T1,'P', P1, self.species)
-                #     dd = self.rho[i] - rho1  #NMOL-nn
-                #     P1 = P1 + dd * 1e5
-                #     if m == itermax:
-                #         raise Exception("Iter max exceeded for rho/P")
-                #     while abs(U_start - U) / U_start > 0.00001 and n < itermax:
-                #         n = n + 1
-                #         U = NMOL * PropsSI('HMOLAR', 'P', P1, 'T', T1, self.species) - self.eta * P1 * self.vol  #Q_inner[i]*tstep
-                #         d = U_start - U
-                #         T1 = T1 + 0.1 * d / U_start * T1
-                #         if n == itermax:
-                #             raise Exception("Iter max exceeded for U/T")
-
-                # self.m_iter[i] = dd
-                # self.n_iter[i] = d
-                # self.U_iter[i] = U / NMOL * self.mass_fluid[i]
+                #print("Iteration: ",i," of ",len(self.P))
+                P1, T1 = self.UDproblem(U_end/ self.mass_fluid[i],self.rho[i],self.P[i-1],self.T_fluid[i-1])
+
                 self.P[i] = P1
                 self.T_fluid[i] = T1
 
@@ -443,7 +420,7 @@ def run(self):
             )
             self.U_mass[i] = PropsSI(
                 "U", "T", self.T_fluid[i], "P", self.P[i], self.species
-            )  # (self.mass_fluid[i] * PropsSI('H', 'P', self.P[i], 'T', self.T_fluid[i], self.species) - self.P[i] * self.vol) / self.mass_fluid[i]
+            )  
             cpcv = PropsSI(
                 "CP0MOLAR", "T", self.T_fluid[i], "P", self.P[i], self.species
             ) / PropsSI("CVMOLAR", "T", self.T_fluid[i], "P", self.P[i], self.species)

input.yml

@@ -8,7 +8,7 @@ vessel:
 initial:
   temperature: 288.0
   pressure: 15000000.
-  fluid: "N2" 
+  fluid: "N2"
 calculation:
   type: "energybalance"
   time_step: 0.05

main.py

@@ -5,6 +5,7 @@
 import yaml
 import sys
 from hyddown import HydDown
+import time
 
 if __name__ == "__main__":
     if len(sys.argv) > 1:
@@ -17,7 +18,10 @@
 
 
     hdown=HydDown(input)
+    start = time.time()
     hdown.run()
+    end = time.time()
+    print('Elapsed time: ',end-start,' sec.')
     hdown.verbose=1
     hdown.plot()
     #hdown.plot("docs/img/N2_filling.png")

requirements.txt

@@ -6,3 +6,4 @@ pytest==6.2.3
 Cerberus==1.3.3
 CoolProp==6.4.1
 PyYAML==5.4.1
+scipy==1.6.0