Merge pull request #640 from dlohmeier/hotfix/vdot_extraction_imcompr…

…essible_real_density

Bugfix: result extraction for incompressible media now based on real density

CHANGELOG.rst

@@ -10,6 +10,7 @@ Change Log
 - [CHANGED] switched from setup.py to pyproject.toml
 - [CHANGED] variable "alpha_w_per_m2k" to "u_w_per_m2k"
 - [FIXED] Pressure plot not working for circ pump
+- [FIXED] volume flow rate for incompressible fluids based on real density, thus in this case results are renamed from "vdot_norm_m3_per_s" to "vdot_m3_per_s"
 - [REMOVED] support for Python 3.8 due to EOL
 
 

src/pandapipes/component_models/component_toolbox.py

@@ -183,19 +183,18 @@ def get_mass_flow_at_nodes(net, node_pit, branch_pit, eg_nodes, comp):
 def standard_branch_wo_internals_result_lookup(net):
     required_results_hyd = [
         ("p_from_bar", "p_from"), ("p_to_bar", "p_to"), ("mdot_to_kg_per_s", "mf_to"),
-        ("mdot_from_kg_per_s", "mf_from"), ("vdot_norm_m3_per_s", "vf"), ("lambda", "lambda"),
-        ("reynolds", "reynolds")
+        ("mdot_from_kg_per_s", "mf_from"), ("lambda", "lambda"), ("reynolds", "reynolds")
     ]
     required_results_ht = [("t_from_k", "temp_from"), ("t_to_k", "temp_to"), ("t_outlet_k", "t_outlet")]
 
     if get_fluid(net).is_gas:
         required_results_hyd.extend([
             ("v_from_m_per_s", "v_gas_from"), ("v_to_m_per_s", "v_gas_to"),
             ("v_mean_m_per_s", "v_gas_mean"), ("normfactor_from", "normfactor_from"),
-            ("normfactor_to", "normfactor_to")
+            ("normfactor_to", "normfactor_to"), ("vdot_norm_m3_per_s", "vf")
         ])
     else:
-        required_results_hyd.extend([("v_mean_m_per_s", "v_mps")])
+        required_results_hyd.extend([("v_mean_m_per_s", "v_mps"), ("vdot_m3_per_s", "vf")])
 
     return required_results_hyd, required_results_ht
 

src/pandapipes/component_models/flow_control_component.py

@@ -129,6 +129,6 @@ def get_result_table(cls, net):
                       "normfactor_to"]
         else:
             output = ["v_mean_m_per_s", "p_from_bar", "p_to_bar", "t_from_k", "t_to_k", "t_outlet_k",
-                      "mdot_from_kg_per_s", "mdot_to_kg_per_s", "vdot_norm_m3_per_s", "reynolds",
+                      "mdot_from_kg_per_s", "mdot_to_kg_per_s", "vdot_m3_per_s", "reynolds",
                       "lambda"]
         return output, True

src/pandapipes/component_models/heat_consumer_component.py

@@ -225,7 +225,7 @@ def get_result_table(cls, net):
                       "normfactor_from", "normfactor_to"]
         else:
             output = ["v_mean_m_per_s", "p_from_bar", "p_to_bar", "t_from_k", "t_to_k", "t_outlet_k", "mdot_from_kg_per_s",
-                      "mdot_to_kg_per_s", "vdot_norm_m3_per_s", "reynolds", "lambda"]
+                      "mdot_to_kg_per_s", "vdot_m3_per_s", "reynolds", "lambda"]
         output += ['deltat_k', 'qext_w']
         return output, True
 

src/pandapipes/component_models/heat_exchanger_component.py

@@ -111,6 +111,6 @@ def get_result_table(cls, net):
                       "normfactor_to"]
         else:
             output = ["v_mean_m_per_s", "p_from_bar", "p_to_bar", "t_from_k", "t_to_k", "t_outlet_k",
-                      "mdot_from_kg_per_s", "mdot_to_kg_per_s", "vdot_norm_m3_per_s", "reynolds",
+                      "mdot_from_kg_per_s", "mdot_to_kg_per_s", "vdot_m3_per_s", "reynolds",
                       "lambda"]
         return output, True

src/pandapipes/component_models/pipe_component.py

@@ -133,18 +133,17 @@ def extract_results(cls, net, options, branch_results, mode):
         res_nodes_from_ht = [("t_from_k", "temp_from")]
         res_nodes_to_hyd = [("p_to_bar", "p_to"), ("mdot_to_kg_per_s", "mf_to")]
         res_nodes_to_ht = [("t_to_k", "temp_to")]
-        res_mean_hyd = [("vdot_norm_m3_per_s", "vf"), ("lambda", "lambda"),
-                        ("reynolds", "reynolds")]
+        res_mean_hyd = [("lambda", "lambda"), ("reynolds", "reynolds")]
         res_branch_ht = [("t_outlet_k", "t_outlet")]
 
         if get_fluid(net).is_gas:
             res_nodes_from_hyd.extend([("v_from_m_per_s", "v_gas_from"),
                                        ("normfactor_from", "normfactor_from")])
             res_nodes_to_hyd.extend([("v_to_m_per_s", "v_gas_to"),
                                      ("normfactor_to", "normfactor_to")])
-            res_mean_hyd.extend([("v_mean_m_per_s", "v_gas_mean")])
+            res_mean_hyd.extend([("v_mean_m_per_s", "v_gas_mean"), ("vdot_norm_m3_per_s", "vf")])
         else:
-            res_mean_hyd.extend([("v_mean_m_per_s", "v_mps")])
+            res_mean_hyd.extend([("v_mean_m_per_s", "v_mps"), ("vdot_m3_per_s", "vf")])
 
         if np.any(cls.get_internal_pipe_number(net) > 1):
             extract_branch_results_with_internals(
@@ -300,7 +299,7 @@ def get_result_table(cls, net):
                       "normfactor_to"]
         else:
             output = ["v_mean_m_per_s", "p_from_bar", "p_to_bar", "t_from_k", "t_to_k", "t_outlet_k",
-                      "mdot_from_kg_per_s", "mdot_to_kg_per_s", "vdot_norm_m3_per_s", "reynolds",
+                      "mdot_from_kg_per_s", "mdot_to_kg_per_s", "vdot_m3_per_s", "reynolds",
                       "lambda"]
         return output, True
 

src/pandapipes/component_models/pressure_control_component.py

@@ -96,17 +96,18 @@ def extract_results(cls, net, options, branch_results, mode):
         """
         required_results_hyd = [
             ("p_from_bar", "p_from"), ("p_to_bar", "p_to"), ("mdot_from_kg_per_s", "mf_from"),
-            ("mdot_to_kg_per_s", "mf_to"), ("vdot_norm_m3_per_s", "vf")
+            ("mdot_to_kg_per_s", "mf_to")
         ]
         required_results_ht = [("t_from_k", "temp_from"), ("t_to_k", "temp_to"), ("t_outlet_k", "t_outlet")]
 
         if get_fluid(net).is_gas:
             required_results_hyd.extend([
                 ("v_from_m_per_s", "v_gas_from"), ("v_to_m_per_s", "v_gas_to"),
-                ("normfactor_from", "normfactor_from"), ("normfactor_to", "normfactor_to")
+                ("normfactor_from", "normfactor_from"), ("normfactor_to", "normfactor_to"),
+                ("vdot_norm_m3_per_s", "vf")
             ])
         else:
-            required_results_hyd.extend([("v_mean_m_per_s", "v_mps")])
+            required_results_hyd.extend([("v_mean_m_per_s", "v_mps"), ("vdot_m3_per_s", "vf")])
 
         extract_branch_results_without_internals(net, branch_results, required_results_hyd,
                                                  required_results_ht, cls.table_name(), mode)
@@ -154,6 +155,6 @@ def get_result_table(cls, net):
                       "vdot_norm_m3_per_s", "normfactor_from", "normfactor_to"]
         else:
             output = ["v_mean_m_per_s", "p_from_bar", "p_to_bar", "t_from_k", "t_to_k", "t_outlet_k",
-                      "mdot_from_kg_per_s", "mdot_to_kg_per_s", "vdot_norm_m3_per_s"]
+                      "mdot_from_kg_per_s", "mdot_to_kg_per_s", "vdot_m3_per_s"]
         output += ["deltap_bar"]
         return output, True

src/pandapipes/component_models/pump_component.py

@@ -138,17 +138,18 @@ def extract_results(cls, net, options, branch_results, mode):
 
         required_results_hyd = [
             ("p_from_bar", "p_from"), ("p_to_bar", "p_to"), ("mdot_to_kg_per_s", "mf_to"),
-            ("mdot_from_kg_per_s", "mf_from"), ("vdot_norm_m3_per_s", "vf"), ("deltap_bar", "pl"),
+            ("mdot_from_kg_per_s", "mf_from"), ("deltap_bar", "pl"),
         ]
         required_results_ht = [("t_from_k", "temp_from"), ("t_to_k", "temp_to"), ("t_outlet_k", "t_outlet")]
 
         if get_fluid(net).is_gas:
             required_results_hyd.extend([
                 ("v_from_m_per_s", "v_gas_from"), ("v_to_m_per_s", "v_gas_to"),
-                ("normfactor_from", "normfactor_from"), ("normfactor_to", "normfactor_to")
+                ("normfactor_from", "normfactor_from"), ("normfactor_to", "normfactor_to"),
+                ("vdot_norm_m3_per_s", "vf")
             ])
         else:
-            required_results_hyd.extend([("v_mean_m_per_s", "v_mps")])
+            required_results_hyd.extend([("v_mean_m_per_s", "v_mps"), ("vdot_m3_per_s", "vf")])
 
         extract_branch_results_without_internals(net, branch_results, required_results_hyd,
                                                  required_results_ht, cls.table_name(), mode)
@@ -224,7 +225,7 @@ def get_result_table(cls, net):
             # TODO: inwieweit sind diese Angaben bei imaginärem Durchmesser sinnvoll?
         else:
             output = ["deltap_bar", "v_mean_m_per_s", "p_from_bar", "p_to_bar", "t_from_k",
-                      "t_to_k", "t_outlet_k", "mdot_from_kg_per_s", "mdot_to_kg_per_s", "vdot_norm_m3_per_s"]
+                      "t_to_k", "t_outlet_k", "mdot_from_kg_per_s", "mdot_to_kg_per_s", "vdot_m3_per_s"]
         if calc_compr_pow:
             output += ["compr_power_mw"]
 

src/pandapipes/component_models/valve_component.py

@@ -91,6 +91,6 @@ def get_result_table(cls, net):
                       "normfactor_to"]
         else:
             output = ["v_mean_m_per_s", "p_from_bar", "p_to_bar", "t_from_k", "t_to_k", "t_outlet_k",
-                      "mdot_from_kg_per_s", "mdot_to_kg_per_s", "vdot_norm_m3_per_s", "reynolds",
+                      "mdot_from_kg_per_s", "mdot_to_kg_per_s", "vdot_m3_per_s", "reynolds",
                       "lambda"]
         return output, True

src/pandapipes/converter/stanet/valve_pipe_component/valve_pipe_component.py

@@ -116,7 +116,7 @@ def get_result_table(cls, net):
                       "p_to_mw",
                       "mdot_from_kg_per_s",
                       "mdot_to_kg_per_s",
-                      "vdot_norm_m3_per_s",
+                      "vdot_m3_per_s",
                       "reynolds",
                       "lambda",
                       "loading_percent"]

src/pandapipes/pf/derivative_calculation.py

@@ -31,7 +31,6 @@ def calculate_derivatives_hydraulic(net, branch_pit, node_pit, options):
     friction_model = options["friction_model"]
     rho = get_branch_real_density(fluid, node_pit, branch_pit)
     eta = get_branch_real_eta(fluid, node_pit, branch_pit)
-    rho_n = fluid.get_density([NORMAL_TEMPERATURE] * len(branch_pit))
 
     lambda_, re = calc_lambda(
         branch_pit[:, MDOTINIT], eta, branch_pit[:, D],
@@ -54,7 +53,7 @@ def calculate_derivatives_hydraulic(net, branch_pit, node_pit, options):
                 as derivatives_hydraulic_incomp
 
         load_vec, load_vec_nodes, df_dm, df_dm_nodes, df_dp, df_dp1 = derivatives_hydraulic_incomp(
-            branch_pit, der_lambda, p_init_i_abs, p_init_i1_abs, height_difference, rho, rho_n)
+            branch_pit, der_lambda, p_init_i_abs, p_init_i1_abs, height_difference, rho)
     else:
         if options["use_numba"]:
             from pandapipes.pf.derivative_toolbox_numba import derivatives_hydraulic_comp_numba \
@@ -65,6 +64,7 @@ def calculate_derivatives_hydraulic(net, branch_pit, node_pit, options):
                 as derivatives_hydraulic_comp, calc_medium_pressure_with_derivative_np as \
                 calc_medium_pressure_with_derivative
         p_m, der_p_m, der_p_m1 = calc_medium_pressure_with_derivative(p_init_i_abs, p_init_i1_abs)
+        rho_n = np.full(len(branch_pit), fluid.get_density(NORMAL_TEMPERATURE))
         comp_fact = fluid.get_compressibility(p_m)
         # TODO: this might not be required
         der_comp = fluid.get_der_compressibility() * der_p_m

src/pandapipes/pf/derivative_toolbox.py

@@ -13,13 +13,16 @@
 
 
 def derivatives_hydraulic_incomp_np(branch_pit, der_lambda, p_init_i_abs, p_init_i1_abs,
-                                    height_difference, rho, rho_n):
+                                    height_difference, rho):
+    # Formulas for pressure loss in incompressible flow
+    # Use medium density ((rho_from + rho_to) / 2) for Darcy Weisbach according to
+    # https://www.schweizer-fn.de/rohr/rohrleitung/rohrleitung.php#fluessigkeiten
     m_init_abs = np.abs(branch_pit[:, MDOTINIT])
     m_init2 = m_init_abs * branch_pit[:, MDOTINIT]
     p_diff = p_init_i_abs - p_init_i1_abs
     const_height = rho * GRAVITATION_CONSTANT * height_difference / P_CONVERSION
     friction_term = np.divide(branch_pit[:, LENGTH] * branch_pit[:, LAMBDA], branch_pit[:, D]) + branch_pit[:, LC]
-    const_term = np.divide(1, branch_pit[:, AREA] ** 2 * rho_n * P_CONVERSION * 2)
+    const_term = np.divide(1, branch_pit[:, AREA] ** 2 * rho * P_CONVERSION * 2)
 
     df_dm = - const_term * (2 * m_init_abs * friction_term + der_lambda
                             * np.divide(branch_pit[:, LENGTH], branch_pit[:, D]) * m_init2)

src/pandapipes/pf/derivative_toolbox_numba.py

@@ -15,9 +15,9 @@
     from numpy import int32, float64, int64
 
 
-@jit((float64[:, :], float64[:], float64[:], float64[:], float64[:], float64[:], float64[:]), nopython=True, cache=False)
+@jit((float64[:, :], float64[:], float64[:], float64[:], float64[:], float64[:]), nopython=True, cache=False)
 def derivatives_hydraulic_incomp_numba(branch_pit, der_lambda, p_init_i_abs, p_init_i1_abs,
-                                       height_difference, rho, rho_n):
+                                       height_difference, rho):
     le = der_lambda.shape[0]
     load_vec = np.zeros_like(der_lambda)
     df_dm = np.zeros_like(der_lambda)
@@ -33,7 +33,7 @@ def derivatives_hydraulic_incomp_numba(branch_pit, der_lambda, p_init_i_abs, p_i
         const_height = rho[i] * GRAVITATION_CONSTANT * height_difference[i] / P_CONVERSION
         friction_term = np.divide(branch_pit[i][LENGTH] * branch_pit[i][LAMBDA], branch_pit[i][D]) \
             + branch_pit[i][LC]
-        const_term = np.divide(1, branch_pit[i][AREA] ** 2 * rho_n[i] * P_CONVERSION * 2)
+        const_term = np.divide(1, branch_pit[i][AREA] ** 2 * rho[i] * P_CONVERSION * 2)
 
         df_dm[i] = -1. * const_term * (2 * m_init_abs * friction_term + der_lambda[i]
                                    * np.divide(branch_pit[i][LENGTH], branch_pit[i][D]) * m_init2)

src/pandapipes/pf/result_extraction.py

@@ -7,6 +7,7 @@
 from pandapipes.pf.internals_toolbox import _sum_by_group
 from pandapipes.pf.pipeflow_setup import get_table_number, get_lookup, get_net_option
 from pandapipes.properties.fluids import get_fluid
+from pandapipes.properties.properties_toolbox import get_branch_real_density
 
 try:
     from numba import jit
@@ -54,11 +55,14 @@ def extract_all_results(net, calculation_mode):
 def get_basic_branch_results(net, branch_pit, node_pit):
     from_nodes = branch_pit[:, FROM_NODE].astype(np.int32)
     to_nodes = branch_pit[:, TO_NODE].astype(np.int32)
-    fluid = get_fluid(net)
     t0 = node_pit[from_nodes, TINIT_NODE]
     t1 = node_pit[to_nodes, TINIT_NODE]
-    vf = branch_pit[:, MDOTINIT] / get_fluid(net).get_density(NORMAL_TEMPERATURE)
-    v = branch_pit[:, MDOTINIT] / fluid.get_density(NORMAL_TEMPERATURE) / branch_pit[:, AREA]
+    fluid = get_fluid(net)
+    if fluid.is_gas:
+        vf = branch_pit[:, MDOTINIT] / fluid.get_density(NORMAL_TEMPERATURE)
+    else:
+        vf = branch_pit[:, MDOTINIT] / get_branch_real_density(fluid, node_pit, branch_pit)
+    v = vf / branch_pit[:, AREA]
     t_outlet = branch_pit[:, TOUTINIT]
     branch_results = {"v_mps": v, "mf_from": branch_pit[:, MDOTINIT], "mf_to": -branch_pit[:, MDOTINIT],
                       "vf": vf, "p_from": node_pit[from_nodes, PINIT], "p_to": node_pit[to_nodes, PINIT],