From acc9b05b52c76beb597e737b71675b120be80366 Mon Sep 17 00:00:00 2001 From: Jolando Kisse Date: Mon, 6 Jan 2025 18:37:44 +0100 Subject: [PATCH 1/5] add documentation for Hofer, add transition area approach --- CHANGELOG.rst | 6 ++- doc/source/components/pipe/pipe_component.rst | 33 +++++++++++--- doc/source/pipeflow/calculation_modes.rst | 2 +- doc/source/references.bib | 38 +++++++++++++++- src/pandapipes/pf/derivative_calculation.py | 43 ++++++++++++++----- src/pandapipes/pf/derivative_toolbox.py | 32 +++++++++++--- src/pandapipes/pf/derivative_toolbox_numba.py | 33 ++++++++++++-- 7 files changed, 157 insertions(+), 30 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 49fc6afcc..52b4374fd 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -1,6 +1,10 @@ Change Log ============= +[upcoming release] - 2025-xx-xx +------------------------------- +- [ADDED] Hofer formula for friction factor "lambda" + [0.11.0] - 2024-11-07 ------------------------------- - [ADDED] heat_consumer plotting @@ -28,8 +32,6 @@ Change Log - [FIXED] alpha also applied to mdot - [REMOVED] support for Python 3.8 due to EOL - - [0.10.0] - 2024-04-09 ------------------------------- diff --git a/doc/source/components/pipe/pipe_component.rst b/doc/source/components/pipe/pipe_component.rst index 95f9e32e2..e3ea6a054 100644 --- a/doc/source/components/pipe/pipe_component.rst +++ b/doc/source/components/pipe/pipe_component.rst @@ -142,13 +142,15 @@ pipe sections. Friction models ^^^^^^^^^^^^^^^ -Three friction models are used to calculate the velocity dependent friction factor: +For friction models are available to calculate the velocity dependent friction factor :math:`\lambda`: -- Nikuradse -- Prandtl-Colebrook -- Swamee-Jain +- Nikuradse ("nikuradse") +- Prandtl-Colebrook ("colebrook") +- Hofer ("hofer") +- Swamee-Jain ("swamee-jain") -Nikuradse is chosen by default. In this case, the friction factor is calculated by: +They are set by the :code:`friction_model` parameter and the name given in parentheses. +*Nikuradse* is chosen by default. In this case, the friction factor is calculated by: .. math:: :nowrap: @@ -161,7 +163,7 @@ Nikuradse is chosen by default. In this case, the friction factor is calculated Note that in literature, Nikuradse is known as a model for turbulent flows. In pandapipes, the formula for the Nikuradse model is also applied for laminar flow. -If Prandtl-Colebrook is selected, the friction factor is calculated iteratively according to +If *Prandtl-Colebrook* (also known as Colebrook-White) is selected, the friction factor is calculated iteratively according to .. math:: :nowrap: @@ -173,7 +175,24 @@ If Prandtl-Colebrook is selected, the friction factor is calculated iteratively Equations for pressure losses due to friction were taken from :cite:`Eberhard1990` and :cite:`Cerbe2008`. -The equation according to Swamee-Jain :cite:`Swamee1976` is an approximation of the calculation method according +The *Hofer-Equation* is an explicit approximation of the Prandtl-Colebrook method and defined as shown in +:cite:`Benner.2019` (based on :cite:`Hofer.1973`) + +.. math:: + :nowrap: + + \begin{align*} + \lambda_H &= \frac{1}{(-2 \log (\frac{4.518}{Re} \cdot \log (\frac{Re}{7}) + \frac{k}{3.71 \cdot d}))^2}\\ + \end{align*} + +Very small Reynolds numbers can lead to problems with the two nested log-functions. +Thus, in pandapipes, the laminar :math:`\lambda_l = 64/Re` is applied for small Reynolds numbers +(:math:`Re < \underline{Re}`) and :math:`\lambda_H` is used for high Reynolds numbers (:math:`Re > \overline{Re}`). +In the transition area :math:`\underline{Re} \leq Re \leq \overline{Re}`, :math:`\lambda_H` is interpolated linearly. +By default, the boundaries of this transition area are set to :math:`\underline{Re}=2000` and :math:`\overline{Re}=3000`. +Note that the derivative used for Hofer is the same as for Prandtl-Colebrook. + +The equation according to *Swamee-Jain* :cite:`Swamee1976` is another approximation of the calculation method according to Prandtl-Colebrook. It is an explicit formula for the friction factor of the transition zone of turbulent flows in pipes and is defined as follows: diff --git a/doc/source/pipeflow/calculation_modes.rst b/doc/source/pipeflow/calculation_modes.rst index b94dc4a81..2994a2f3a 100644 --- a/doc/source/pipeflow/calculation_modes.rst +++ b/doc/source/pipeflow/calculation_modes.rst @@ -92,7 +92,7 @@ In gas flows, the velocity is typically not constant along a pipeline. For this tables for pipes show more entries in comparison with the result tables for incompressible media. -Temperature calculations (pipeflow option: mode = "sequential", mode = "bidrectional" or mode = "heat") +Temperature calculations (pipeflow option: mode = "sequential", mode = "bidirectional" or mode = "heat") ======================================================================================================= Important parameters of the network main components (junctions and pipes) needed for the calculation diff --git a/doc/source/references.bib b/doc/source/references.bib index a09cda44d..500191f78 100644 --- a/doc/source/references.bib +++ b/doc/source/references.bib @@ -81,4 +81,40 @@ @article{Schmidt.2015 issn = {1389-4420}, journal = {Optimization and Engineering}, doi = {10.1007/s11081-014-9246-x} -} \ No newline at end of file +} + + +@incollection{Benner.2019, + author = {Benner, Peter and Grundel, Sara and Himpe, Christian and Huck, Christoph and Streubel, Tom and Tischendorf, Caren}, + title = {Gas Network Benchmark Models}, + pages = {171--197}, + publisher = {{Springer International Publishing}}, + isbn = {978-3-030-03717-8}, + series = {Differential-Algebraic Equations Forum}, + editor = {Campbell, Stephen and Ilchmann, Achim and Mehrmann, Volker and Reis, Timo}, + booktitle = {Applications of Differential-Algebraic Equations: Examples and Benchmarks}, + year = {2019}, + address = {Cham}, + doi = {10.1007/11221{\_}2018{\_}5}, +} + +@book{Campbell.2019, + year = {2019}, + title = {Applications of Differential-Algebraic Equations: Examples and Benchmarks}, + address = {Cham}, + publisher = {{Springer International Publishing}}, + isbn = {978-3-030-03717-8}, + series = {Differential-Algebraic Equations Forum}, + editor = {Campbell, Stephen and Ilchmann, Achim and Mehrmann, Volker and Reis, Timo}, + doi = {10.1007/978-3-030-03718-5} +} + +@article{Hofer.1973, + author = {Hofer, P.}, + year = {1973}, + title = {Beurteilung von Fehlern in Rohrnetzberechnungen}, + pages = {113--119}, + volume = {114}, + number = {3}, + journal = {Das Gas- und Wasserfach. GWF - Ausgabe Wasser, Abwasser} +} diff --git a/src/pandapipes/pf/derivative_calculation.py b/src/pandapipes/pf/derivative_calculation.py index 922f5b09b..891287ab5 100644 --- a/src/pandapipes/pf/derivative_calculation.py +++ b/src/pandapipes/pf/derivative_calculation.py @@ -140,6 +140,20 @@ def get_derived_values(node_pit, from_nodes, to_nodes, use_numba): return calc_derived_values_np(node_pit, from_nodes, to_nodes) +def calc_lambda_transition_area(lambda_laminar, lambda_turb, re, begin_transition_re=2000, + end_transition_re=3000): + """ + Linear interpolation in the transition area between laminar and turbulent flow. Critical + point is usually at Reynolds = 2300, thus the default transition area is around that point. + """ + lambda_tot = lambda_laminar + share_turb = (re - begin_transition_re) / (end_transition_re - begin_transition_re) + lambda_transition = (1 - share_turb) * lambda_laminar + share_turb * lambda_turb + lambda_tot[(re >= begin_transition_re) & (re < end_transition_re)] = \ + lambda_transition[(re >= begin_transition_re) & (re < end_transition_re)] + lambda_tot[re >= end_transition_re] = lambda_turb[re >= end_transition_re] + return lambda_tot + def calc_lambda(m, eta, d, k, gas_mode, friction_model, lengths, options, area): """ Function calculates the friction factor of a pipe. Turbulence is calculated based on @@ -171,16 +185,19 @@ def calc_lambda(m, eta, d, k, gas_mode, friction_model, lengths, options, area): from pandapipes.pf.derivative_toolbox_numba import calc_lambda_nikuradse_incomp_numba as \ calc_lambda_nikuradse_incomp, colebrook_numba as colebrook, \ calc_lambda_nikuradse_comp_numba as calc_lambda_nikuradse_comp + from pandapipes.pf.derivative_toolbox import calc_lambda_hofer_comp_np as calc_lambda_hofer_comp # numba version + # not yet implemented else: from pandapipes.pf.derivative_toolbox import calc_lambda_nikuradse_incomp_np as \ calc_lambda_nikuradse_incomp, colebrook_np as colebrook, \ - calc_lambda_nikuradse_comp_np as calc_lambda_nikuradse_comp + calc_lambda_nikuradse_comp_np as calc_lambda_nikuradse_comp, \ + calc_lambda_hofer_comp_np as calc_lambda_hofer_comp if gas_mode: re, lambda_laminar, lambda_nikuradse = calc_lambda_nikuradse_comp(m, d, k, eta, area) else: re, lambda_laminar, lambda_nikuradse = calc_lambda_nikuradse_incomp(m, d, k, eta, area) - if friction_model == "colebrook": + if friction_model.lower() == "colebrook": # TODO: move this import to top level if possible from pandapipes.pipeflow import PipeflowNotConverged max_iter = options.get("max_iter_colebrook", 100) @@ -192,9 +209,14 @@ def calc_lambda(m, eta, d, k, gas_mode, friction_model, lengths, options, area): "inconsistencies. The maximum iterations can be given as 'max_iter_colebrook' " "argument to the pipeflow.") return lambda_colebrook, re - elif friction_model == "swamee-jain": + elif friction_model.lower() == "swamee-jain": lambda_swamee_jain = 0.25 / ((np.log10(k / (3.7 * d) + 5.74 / (re ** 0.9))) ** 2) return lambda_swamee_jain, re + elif friction_model.lower() == "hofer": + re, lambda_hofer = calc_lambda_hofer_comp(m, d, k, eta, area) + lambda_tot = calc_lambda_transition_area(lambda_laminar, lambda_hofer, re, + begin_transition_re=2000, end_transition_re=4000) + return lambda_tot, re else: # lambda_tot = np.where(re > 2300, lambda_laminar + lambda_nikuradse, lambda_laminar) lambda_tot = lambda_laminar + lambda_nikuradse @@ -203,12 +225,12 @@ def calc_lambda(m, eta, d, k, gas_mode, friction_model, lengths, options, area): def calc_der_lambda(m, eta, d, k, friction_model, lambda_pipe, area): """ - Function calculates the derivative of lambda with respect to v. Turbulence is calculated based - on Nikuradse. This should not be a problem as the pressure loss term will equal zero - (lambda * u^2). + Function calculates the derivative of lambda with respect to m (mass flow rate). Turbulence is + calculated based on Nikuradse. This should not be a problem as the pressure loss term will + equal zero (lambda * u^2). - :param v: - :type v: + :param m: + :type m: :param eta: :type eta: :param rho: @@ -231,7 +253,8 @@ def calc_der_lambda(m, eta, d, k, friction_model, lambda_pipe, area): lambda_der = np.zeros_like(m) pos = m != 0 - if friction_model == "colebrook": + if friction_model.lower() in ["colebrook", "hofer"]: # hofer is explicit approximation of + # colebrook b_term[pos] = (2.51 * eta[pos] * area[pos] / (m[pos] * d[pos] * np.sqrt(lambda_pipe[pos])) + k[pos] / (3.71 * d[pos])) @@ -244,7 +267,7 @@ def calc_der_lambda(m, eta, d, k, friction_model, lambda_pipe, area): lambda_der[pos] = df_dm[pos] / df_dlambda[pos] return lambda_der - elif friction_model == "swamee-jain": + elif friction_model.lower() == "swamee-jain": param = (k[pos] / (3.7 * d[pos]) + 5.74 * ((eta[pos] * area[pos]) / (np.abs(m[pos]) * d[pos])) ** 0.9) # 0.5 / (log(10) * log(param)^3 * param) * 5.166 * abs(eta)^0.9 / (abs(rho * d)^0.9 diff --git a/src/pandapipes/pf/derivative_toolbox.py b/src/pandapipes/pf/derivative_toolbox.py index d627560f3..330a52495 100644 --- a/src/pandapipes/pf/derivative_toolbox.py +++ b/src/pandapipes/pf/derivative_toolbox.py @@ -4,6 +4,7 @@ import numpy as np from numpy import linalg +from copy import deepcopy from pandapipes.constants import P_CONVERSION, GRAVITATION_CONSTANT, NORMAL_PRESSURE, \ NORMAL_TEMPERATURE @@ -73,24 +74,45 @@ def derivatives_hydraulic_comp_np(node_pit, branch_pit, lambda_, der_lambda, p_i return load_vec, load_vec_nodes_from, load_vec_nodes_to, df_dm, df_dm_nodes, df_dp, df_dp1 -def calc_lambda_nikuradse_incomp_np(m, d, k, eta, area): +def calc_lambda_laminar(m, d, eta, area): m_abs = np.abs(m) re = np.divide(m_abs * d, eta * area) lambda_laminar = np.zeros_like(m) lambda_laminar[~np.isclose(re, 0)] = 64 / re[~np.isclose(re, 0)] + return lambda_laminar, re + + +def calc_lambda_nikuradse_incomp_np(m, d, k, eta, area): + lambda_laminar, re = calc_lambda_laminar(m, d, eta, area) lambda_nikuradse = np.divide(1, (-2 * np.log10(k / (3.71 * d))) ** 2) return re, lambda_laminar, lambda_nikuradse def calc_lambda_nikuradse_comp_np(m, d, k, eta, area): - m_abs = np.abs(m) - re = np.divide(m_abs * d, eta * area) - lambda_laminar = np.zeros_like(m) - lambda_laminar[~np.isclose(re, 0)] = 64 / re[~np.isclose(re, 0)] + lambda_laminar, re = calc_lambda_laminar(m, d, eta, area) lambda_nikuradse = np.divide(1, (2 * np.log10(d / k) + 1.14) ** 2) return re, lambda_laminar, lambda_nikuradse +def calc_lambda_hofer_comp_np(m, d, k, eta, area, hofer_re_threshold=2000): + """Calculate Lambda acc. to Hofer-Equation (explicit version of Colebrook) [1]. + Due to nested log10 operations, small Reynolds numbers must be avoided. + It is justified because the laminar Lambda (instead of Hofer) will be applied for small + Reynolds numbers anyway. + + [1]: P. Hofer. Beurteilung von Fehlern in Rohrnetzberechnungen (Error evaluation in calculation + of pipelines). GWF–Gas/Erdgas, 114(3):113–119, 1973 + """ + lambda_laminar, re = calc_lambda_laminar(m, d, eta, area) + lambda_hofer = np.zeros_like(re) + re_lower_lim = deepcopy(re) + re_lower_lim[re < hofer_re_threshold] = hofer_re_threshold + lambda_hofer = \ + np.divide(1, (-2 * np.log10((4.518/re_lower_lim) * np.log10(re_lower_lim/7) + (k / (3.71 * d)))) ** 2) + lambda_hofer[re < hofer_re_threshold] = 0 + return re, lambda_hofer + + def calc_medium_pressure_with_derivative_np(p_init_i_abs, p_init_i1_abs): val = 2 / 3 p_m = p_init_i_abs.copy() diff --git a/src/pandapipes/pf/derivative_toolbox_numba.py b/src/pandapipes/pf/derivative_toolbox_numba.py index 69be5c7fd..9ff02d4bd 100644 --- a/src/pandapipes/pf/derivative_toolbox_numba.py +++ b/src/pandapipes/pf/derivative_toolbox_numba.py @@ -111,15 +111,40 @@ def calc_lambda_nikuradse_comp_numba(m, d, k, eta, area): lambda_nikuradse = np.empty_like(m) lambda_laminar = np.zeros_like(m) re = np.empty_like(m) + m_abs = np.abs(m) for i, mi in enumerate(m): - m_abs = np.abs(mi) - re[i] = np.divide(m_abs * d[i], eta[i] * area[i]) + re[i] = (m_abs[i] * d[i]) / (eta[i] * area[i]) if re[i] != 0: - lambda_laminar[i] = np.divide(64, re[i]) - lambda_nikuradse[i] = np.divide(1, (2 * np.log10(np.divide(d[i], k[i])) + 1.14) ** 2) + lambda_laminar[i] = 64 / re[i] + lambda_nikuradse[i] = (2 * np.log10((d[i] / k[i])) + 1.14) ** -2 return re, lambda_laminar, lambda_nikuradse +@jit((float64[:], float64[:], float64[:], float64[:], float64[:], int64), nopython=True) +def calc_lambda_hofer_comp_numba(m, d, k, eta, area, hofer_re_threshold=2000): + """Calculate Lambda acc. to Hofer-Equation (explicit version of Colebrook) [1]. + Due to nested log10 operations, small Reynolds numbers must be avoided. + It is justified because the laminar Lambda (instead of Hofer) will be applied for small + Reynolds numbers anyway. + + [1]: P. Hofer. Beurteilung von Fehlern in Rohrnetzberechnungen (Error evaluation in calculation + of pipelines). GWF–Gas/Erdgas, 114(3):113–119, 1973 + """ + lambda_hofer = np.zeros_like(m) + lambda_laminar = np.zeros_like(m) + re = np.empty_like(m) + m_abs = np.abs(m) + for i, mi in enumerate(m): + re[i] = (m_abs[i] * d[i]) / (eta[i] * area[i]) + if re[i] != 0: + lambda_laminar[i] = 64 / re[i] + if re[i] >= hofer_re_threshold: + lambda_hofer[i] = (-2 * np.log10((4.518/re[i]) * np.log10(re[i]/7) + + (k[i] / (3.71 * d[i]))) + ) ** -2 + return re, lambda_laminar, lambda_hofer + + @jit((float64[:], float64[:]), nopython=True, cache=False) def calc_medium_pressure_with_derivative_numba(p_init_i_abs, p_init_i1_abs): p_m = p_init_i_abs.copy() From 6fae5f6002a5a92858349b37447fd5955365389b Mon Sep 17 00:00:00 2001 From: Jolando Kisse Date: Mon, 6 Jan 2025 18:49:04 +0100 Subject: [PATCH 2/5] improve doc: dynamic parantheses size --- doc/source/components/pipe/pipe_component.rst | 8 ++++---- doc/source/references.bib | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/doc/source/components/pipe/pipe_component.rst b/doc/source/components/pipe/pipe_component.rst index e3ea6a054..2115be8b9 100644 --- a/doc/source/components/pipe/pipe_component.rst +++ b/doc/source/components/pipe/pipe_component.rst @@ -156,7 +156,7 @@ They are set by the :code:`friction_model` parameter and the name given in paren :nowrap: \begin{align*} - \lambda &= \frac{64}{Re} + \frac{1}{(-2 \cdot \log (\frac{k}{3.71 \cdot d}))^2}\\ + \lambda &= \frac{64}{Re} + \frac{1}{\left(-2 \cdot \log \left(\frac{k}{3.71 \cdot d}\right)\right)^2}\\ \end{align*} @@ -169,7 +169,7 @@ If *Prandtl-Colebrook* (also known as Colebrook-White) is selected, the friction :nowrap: \begin{align*} - \frac{1}{\sqrt{\lambda}} &= -2 \cdot \log (\frac{2.51}{Re \cdot \sqrt{\lambda}} + \frac{k}{3.71 \cdot d})\\ + \frac{1}{\sqrt{\lambda}} &= -2 \cdot \log \left(\frac{2.51}{Re \cdot \sqrt{\lambda}} + \frac{k}{3.71 \cdot d}\right)\\ \end{align*} Equations for pressure losses due to friction were taken from :cite:`Eberhard1990` and @@ -182,7 +182,7 @@ The *Hofer-Equation* is an explicit approximation of the Prandtl-Colebrook metho :nowrap: \begin{align*} - \lambda_H &= \frac{1}{(-2 \log (\frac{4.518}{Re} \cdot \log (\frac{Re}{7}) + \frac{k}{3.71 \cdot d}))^2}\\ + \lambda_H &= \frac{1}{(-2 \log \left(\frac{4.518}{Re} \cdot \log \left(\frac{Re}{7}) + \frac{k}{3.71 \cdot d}\right)\right)^2}\\ \end{align*} Very small Reynolds numbers can lead to problems with the two nested log-functions. @@ -200,7 +200,7 @@ zone of turbulent flows in pipes and is defined as follows: :nowrap: \begin{align*} - \lambda &= \frac{0.25}{(\log(\frac{k}{3.7 \cdot d} + \frac{5.74}{Re^{0.9}}))^2}\\ + \lambda &= \frac{0.25}{\left(\log\left(\frac{k}{3.7 \cdot d} + \frac{5.74}{Re^{0.9}}\right)\right)^2}\\ \end{align*} diff --git a/doc/source/references.bib b/doc/source/references.bib index 500191f78..9dbb1f6d3 100644 --- a/doc/source/references.bib +++ b/doc/source/references.bib @@ -95,7 +95,7 @@ @incollection{Benner.2019 booktitle = {Applications of Differential-Algebraic Equations: Examples and Benchmarks}, year = {2019}, address = {Cham}, - doi = {10.1007/11221{\_}2018{\_}5}, + doi = {10.1007/11221_2018_5}, } @book{Campbell.2019, @@ -112,7 +112,7 @@ @book{Campbell.2019 @article{Hofer.1973, author = {Hofer, P.}, year = {1973}, - title = {Beurteilung von Fehlern in Rohrnetzberechnungen}, + title = {{Beurteilung von Fehlern in Rohrnetzberechnungen}}, pages = {113--119}, volume = {114}, number = {3}, From 0ce86f89c3129ced72129374190f9cc4a3dde59d Mon Sep 17 00:00:00 2001 From: Jolando Kisse Date: Tue, 7 Jan 2025 10:28:41 +0100 Subject: [PATCH 3/5] extend headline underlining --- doc/source/pipeflow/calculation_modes.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/source/pipeflow/calculation_modes.rst b/doc/source/pipeflow/calculation_modes.rst index 2994a2f3a..fd47f60d7 100644 --- a/doc/source/pipeflow/calculation_modes.rst +++ b/doc/source/pipeflow/calculation_modes.rst @@ -93,7 +93,7 @@ tables for pipes show more entries in comparison with the result tables for inco Temperature calculations (pipeflow option: mode = "sequential", mode = "bidirectional" or mode = "heat") -======================================================================================================= +======================================================================================================== Important parameters of the network main components (junctions and pipes) needed for the calculation are listed in the following table. The :ref:`component section ` of this manual contains From c83b746e114b6693789e88e2c094a52e2e765937 Mon Sep 17 00:00:00 2001 From: Jolando Kisse Date: Tue, 7 Jan 2025 16:40:15 +0100 Subject: [PATCH 4/5] improve handling of mdot=0, use correct hofer_numba --- src/pandapipes/pf/derivative_calculation.py | 11 +++++------ src/pandapipes/pf/derivative_toolbox.py | 2 +- src/pandapipes/pf/derivative_toolbox_numba.py | 5 +++-- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/pandapipes/pf/derivative_calculation.py b/src/pandapipes/pf/derivative_calculation.py index 891287ab5..ee12eb123 100644 --- a/src/pandapipes/pf/derivative_calculation.py +++ b/src/pandapipes/pf/derivative_calculation.py @@ -157,7 +157,7 @@ def calc_lambda_transition_area(lambda_laminar, lambda_turb, re, begin_transitio def calc_lambda(m, eta, d, k, gas_mode, friction_model, lengths, options, area): """ Function calculates the friction factor of a pipe. Turbulence is calculated based on - Nikuradse. If v equals 0, a value of 0.001 is used in order to avoid division by zero. + Nikuradse. If m equals 0, a value of 0.001 is used in order to avoid division by zero. This should not be a problem as the pressure loss term will equal zero (lambda * u^2). :param m: @@ -184,9 +184,8 @@ def calc_lambda(m, eta, d, k, gas_mode, friction_model, lengths, options, area): if options["use_numba"]: from pandapipes.pf.derivative_toolbox_numba import calc_lambda_nikuradse_incomp_numba as \ calc_lambda_nikuradse_incomp, colebrook_numba as colebrook, \ - calc_lambda_nikuradse_comp_numba as calc_lambda_nikuradse_comp - from pandapipes.pf.derivative_toolbox import calc_lambda_hofer_comp_np as calc_lambda_hofer_comp # numba version - # not yet implemented + calc_lambda_nikuradse_comp_numba as calc_lambda_nikuradse_comp, \ + calc_lambda_hofer_comp_numba as calc_lambda_hofer_comp else: from pandapipes.pf.derivative_toolbox import calc_lambda_nikuradse_incomp_np as \ calc_lambda_nikuradse_incomp, colebrook_np as colebrook, \ @@ -213,7 +212,7 @@ def calc_lambda(m, eta, d, k, gas_mode, friction_model, lengths, options, area): lambda_swamee_jain = 0.25 / ((np.log10(k / (3.7 * d) + 5.74 / (re ** 0.9))) ** 2) return lambda_swamee_jain, re elif friction_model.lower() == "hofer": - re, lambda_hofer = calc_lambda_hofer_comp(m, d, k, eta, area) + re, lambda_laminar, lambda_hofer = calc_lambda_hofer_comp(m, d, k, eta, area) lambda_tot = calc_lambda_transition_area(lambda_laminar, lambda_hofer, re, begin_transition_re=2000, end_transition_re=4000) return lambda_tot, re @@ -251,7 +250,7 @@ def calc_der_lambda(m, eta, d, k, friction_model, lambda_pipe, area): df_dm = np.zeros_like(m) df_dlambda = np.zeros_like(m) lambda_der = np.zeros_like(m) - pos = m != 0 + pos = ~np.isclose(abs(m), 0) if friction_model.lower() in ["colebrook", "hofer"]: # hofer is explicit approximation of # colebrook diff --git a/src/pandapipes/pf/derivative_toolbox.py b/src/pandapipes/pf/derivative_toolbox.py index 330a52495..ab67c5dfe 100644 --- a/src/pandapipes/pf/derivative_toolbox.py +++ b/src/pandapipes/pf/derivative_toolbox.py @@ -110,7 +110,7 @@ def calc_lambda_hofer_comp_np(m, d, k, eta, area, hofer_re_threshold=2000): lambda_hofer = \ np.divide(1, (-2 * np.log10((4.518/re_lower_lim) * np.log10(re_lower_lim/7) + (k / (3.71 * d)))) ** 2) lambda_hofer[re < hofer_re_threshold] = 0 - return re, lambda_hofer + return re, lambda_laminar, lambda_hofer def calc_medium_pressure_with_derivative_np(p_init_i_abs, p_init_i1_abs): diff --git a/src/pandapipes/pf/derivative_toolbox_numba.py b/src/pandapipes/pf/derivative_toolbox_numba.py index 9ff02d4bd..e91fb0eef 100644 --- a/src/pandapipes/pf/derivative_toolbox_numba.py +++ b/src/pandapipes/pf/derivative_toolbox_numba.py @@ -120,8 +120,8 @@ def calc_lambda_nikuradse_comp_numba(m, d, k, eta, area): return re, lambda_laminar, lambda_nikuradse -@jit((float64[:], float64[:], float64[:], float64[:], float64[:], int64), nopython=True) -def calc_lambda_hofer_comp_numba(m, d, k, eta, area, hofer_re_threshold=2000): +@jit((float64[:], float64[:], float64[:], float64[:], float64[:]), nopython=True) +def calc_lambda_hofer_comp_numba(m, d, k, eta, area): """Calculate Lambda acc. to Hofer-Equation (explicit version of Colebrook) [1]. Due to nested log10 operations, small Reynolds numbers must be avoided. It is justified because the laminar Lambda (instead of Hofer) will be applied for small @@ -130,6 +130,7 @@ def calc_lambda_hofer_comp_numba(m, d, k, eta, area, hofer_re_threshold=2000): [1]: P. Hofer. Beurteilung von Fehlern in Rohrnetzberechnungen (Error evaluation in calculation of pipelines). GWF–Gas/Erdgas, 114(3):113–119, 1973 """ + hofer_re_threshold = 2000 lambda_hofer = np.zeros_like(m) lambda_laminar = np.zeros_like(m) re = np.empty_like(m) From ea07daea2f869dd1dca3a8e72def96565f7085ea Mon Sep 17 00:00:00 2001 From: Jolando Kisse Date: Wed, 15 Jan 2025 14:32:44 +0100 Subject: [PATCH 5/5] further documentation, add deprecation warning for max_iter, increase default max_iter_hyd --- CHANGELOG.rst | 3 +++ doc/source/components/pipe/pipe_component.rst | 2 +- src/pandapipes/pf/derivative_calculation.py | 19 +++++++++---------- src/pandapipes/pf/pipeflow_setup.py | 7 ++++--- src/pandapipes/pipeflow.py | 8 ++++++++ 5 files changed, 25 insertions(+), 14 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 52b4374fd..2bb82f186 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -4,6 +4,9 @@ Change Log [upcoming release] - 2025-xx-xx ------------------------------- - [ADDED] Hofer formula for friction factor "lambda" +- [ADDED] deprecation warning for "max_iter" (use "max_iter_hyd" or "max_iter_therm" instead) +- [CHANGED] default maximum iterations for hydraulic calculations now 30 instead of 10 + [0.11.0] - 2024-11-07 ------------------------------- diff --git a/doc/source/components/pipe/pipe_component.rst b/doc/source/components/pipe/pipe_component.rst index 2115be8b9..8e4c783b5 100644 --- a/doc/source/components/pipe/pipe_component.rst +++ b/doc/source/components/pipe/pipe_component.rst @@ -142,7 +142,7 @@ pipe sections. Friction models ^^^^^^^^^^^^^^^ -For friction models are available to calculate the velocity dependent friction factor :math:`\lambda`: +Four friction models are available to calculate the velocity dependent friction factor :math:`\lambda`: - Nikuradse ("nikuradse") - Prandtl-Colebrook ("colebrook") diff --git a/src/pandapipes/pf/derivative_calculation.py b/src/pandapipes/pf/derivative_calculation.py index ee12eb123..0f4273d38 100644 --- a/src/pandapipes/pf/derivative_calculation.py +++ b/src/pandapipes/pf/derivative_calculation.py @@ -156,24 +156,24 @@ def calc_lambda_transition_area(lambda_laminar, lambda_turb, re, begin_transitio def calc_lambda(m, eta, d, k, gas_mode, friction_model, lengths, options, area): """ - Function calculates the friction factor of a pipe. Turbulence is calculated based on - Nikuradse. If m equals 0, a value of 0.001 is used in order to avoid division by zero. - This should not be a problem as the pressure loss term will equal zero (lambda * u^2). + Function calculates the friction factor of a pipe, considering laminar and turbulent flow. + Turbulence is calculated based on the friction model provided. If m equals 0, a value of 0.001 is used in order + to avoid division by zero. This should not be a problem as the pressure loss term will equal zero (lambda * u^2). - :param m: + :param m: mass flow rate :type m: :param eta: :type eta: :param rho: :type rho: - :param d: + :param d: inner pipe diameter :type d: - :param k: + :param k: roughness coefficient :type k: :param gas_mode: :type gas_mode: - :param friction_model: - :type friction_model: + :param friction_model: friction formula, one of "nikuradse", "colebrook", "hofer", "swamee-jain" + :type friction_model: str :param lengths: :type lengths: :param options: @@ -225,8 +225,7 @@ def calc_lambda(m, eta, d, k, gas_mode, friction_model, lengths, options, area): def calc_der_lambda(m, eta, d, k, friction_model, lambda_pipe, area): """ Function calculates the derivative of lambda with respect to m (mass flow rate). Turbulence is - calculated based on Nikuradse. This should not be a problem as the pressure loss term will - equal zero (lambda * u^2). + calculated based on the friction model provided. :param m: :type m: diff --git a/src/pandapipes/pf/pipeflow_setup.py b/src/pandapipes/pf/pipeflow_setup.py index a905fba59..c34fd1ffc 100644 --- a/src/pandapipes/pf/pipeflow_setup.py +++ b/src/pandapipes/pf/pipeflow_setup.py @@ -33,8 +33,9 @@ logger = logging.getLogger(__name__) -default_options = {"friction_model": "nikuradse", "tol_p": 1e-5, "tol_m": 1e-5, - "tol_T": 1e-3, "tol_res": 1e-3, "max_iter_hyd": 10, "max_iter_therm": 10, "max_iter_bidirect": 10, +default_options = {"friction_model": "nikuradse", + "tol_p": 1e-5, "tol_m": 1e-5, "tol_T": 1e-3, "tol_res": 1e-3, + "max_iter_hyd": 30, "max_iter_therm": 10, "max_iter_bidirect": 10, "error_flag": False, "alpha": 1, "nonlinear_method": "constant", "mode": "hydraulics", "ambient_temperature": 293.15, "check_connectivity": True, @@ -235,7 +236,7 @@ def init_options(net, local_parameters): calculation of the barometric formula - **friction_model** (str): "nikuradse" - The friction model that shall be used to identify\ - the value for lambda (can be "nikuradse" or "colebrook") + the value for lambda (can be "nikuradse", "colebrook", "hofer", or "swamee-jain") - **alpha** (float): 1 - The step width for the Newton iterations. If the Newton steps \ shall be damped, **alpha** can be reduced. See also the **nonlinear_method** \ diff --git a/src/pandapipes/pipeflow.py b/src/pandapipes/pipeflow.py index 4ab146848..4a8673fcd 100644 --- a/src/pandapipes/pipeflow.py +++ b/src/pandapipes/pipeflow.py @@ -4,6 +4,7 @@ import numpy as np from numpy import linalg +import warnings from scipy.sparse.linalg import spsolve from pandapipes.idx_branch import MDOTINIT, TOUTINIT, FROM_NODE_T_SWITCHED @@ -170,6 +171,13 @@ def hydraulics(net): # --------------------------------------------------------------------------------------------- net.converged = False reduce_pit(net, mode="hydraulics") + if hasattr(net, "_options") and "max_iter" in net._options.keys(): + warnings.warn("The net option 'max_iter' is deprecated. Please use 'max_iter_hyd', 'max_iter_therm', " + "or 'max_iter_bidirect' instead.", + FutureWarning) + net._options["max_iter_hyd"] = max(net._options["max_iter_hyd"], net._options["max_iter"]) + net._options["max_iter_therm"] = max(net._options["max_iter_therm"], net._options["max_iter"]) + _ = net._options.pop("max_iter") if not get_net_option(net, "reuse_internal_data") or "_internal_data" not in net: net["_internal_data"] = dict() solver_vars = ['mdot', 'p', 'mdotslack']