Hydraulics of hot water in pipe: Unexpected effects of temperatures on pressure drop

[j-zipplies]

Describe the bug
When simulating hot water in a pipe via pp.pipeflow(net, mode="all"), the results for velocity and re-Number => lambda and in final consequence the pressure drop have unexpected (in)dependencies on the temperatures. Probably, this comes from the way density and viscosity are calculated and the fact that hydraulic and thermal calculation are fully separated.

1. Density
   Apparently the density does not depend at all on the temperatures in the net. It is apparently set at t=273.15 K (999.9 kg/m³).
   In consequence, the flow velocity in res_pipes does not at all depend on the temperatures in the net, it just depends on mdot.

2. Viscosity
   The viscosity seems to be calculated according to tfluid_k at junctions - or, if applicable, the temperature of external grid.
   In consequence, reynolds numbers => lambda => pressure drop in the pipes depend on these temperatures.

To Reproduce
I set up an example, external grid => pipe1 => pipe2 => sink, with four different configurations of temperatures (tfluid_k at junctions and the t_k from the external grid are either 300 or 350).

import math
import pandas as pd
import pandapipes as pp

def example_hot_water():
    net = pp.create_empty_network("network", fluid="water")
    mdot_kg_per_s = 10
    temps_k = [300, 350]
    d_pipe = 0.1
    l_pipe = 1
    
    for t_fluid in temps_k:
        for t_j_k in temps_k:
            create_source_2pipes_sink(net,t_source_k=t_fluid, t_junctions_k=t_j_k,
                                    pipe_name_suffix=f"_tj_{t_j_k}_tf_{t_fluid}",
                                    mdot_kg_per_s=mdot_kg_per_s, diameter_m=d_pipe, length_km=l_pipe)

    pp.pipeflow(net, mode="all", friction_model="colebrook")

    results = pd.DataFrame.join(
        net.pipe.loc[:,["name"]],
        net.res_pipe.loc[:,["t_from_k","v_mean_m_per_s","reynolds","lambda"]]
    )
    results.loc[:,'dp'] = net.res_pipe.p_from_bar - net.res_pipe.p_to_bar

    # Density (for velocity calculation) is calculated at 273.15 K, which is not the same as the actual temperature of the fluid
    velocities_equal = net.res_pipe['v_mean_m_per_s'].nunique() == 1
    print(f"All flow velocities equal: {velocities_equal}")
    rho_from_v_mean = round(mdot_kg_per_s / (net.res_pipe.loc[:,"v_mean_m_per_s"] * math.pi * d_pipe**2 / 4),2)
    print(f"Calculated densities from v_mean: {set(rho_from_v_mean)}")
    for t_k in [273.15]+temps_k:
        print(f"Actual density at {t_k} K: {net.fluid.get_density(t_k)}")
        
    # dp is really calculated from v_mean, the default density and the various lambdas
    dp_calc = net.res_pipe.loc[:,'lambda'] * l_pipe*1000 / d_pipe * 0.5 * rho_from_v_mean * net.res_pipe.loc[:,"v_mean_m_per_s"]**2 /10**5
    results.loc[:,'dp_calc'] = dp_calc
    print(results)

    return net

def create_source_2pipes_sink(net,t_source_k, t_junctions_k, mdot_kg_per_s,  length_km, diameter_m, pipe_name_suffix='', k_mm=0.1):

    j1, j2, j3 = pp.create_junctions(net, nr_junctions=3, pn_bar=1, tfluid_k=t_junctions_k)
    pp.create_pipe_from_parameters(net, from_junction=j1, to_junction=j2, name="pipe1"+pipe_name_suffix,
        length_km=length_km, diameter_m=diameter_m, k_mm=k_mm)
    pp.create_pipe_from_parameters(net, from_junction=j2, to_junction=j3, name="pipe2"+pipe_name_suffix,
        length_km=length_km, diameter_m=diameter_m, k_mm=k_mm)
    pp.create_ext_grid(net, junction=j1, p_bar=10, t_k=t_source_k)
    pp.create_sink(net, junction=j3, mdot_kg_per_s=mdot_kg_per_s)

if __name__ == "__main__":
    net = example_hot_water()

Output:

All flow velocities equal: True
Calculated densities from v_mean: {999.88}
Actual density at 273.15 K: 999.8801666666666
Actual density at 300 K: 996.49
Actual density at 350 K: 973.62
                  name  t_from_k  v_mean_m_per_s       reynolds    lambda  \
0  pipe1_tj_300_tf_300     300.0        1.273392  149231.076504  0.021434   
1  pipe2_tj_300_tf_300     300.0        1.273392  149231.076504  0.021434   
2  pipe1_tj_350_tf_300     300.0        1.273392  240324.564880  0.020810   
3  pipe2_tj_350_tf_300     300.0        1.273392  345613.340048  0.020473   
4  pipe1_tj_300_tf_350     350.0        1.273392  240324.564880  0.020810   
5  pipe2_tj_300_tf_350     350.0        1.273392  149231.076504  0.021434   
6  pipe1_tj_350_tf_350     350.0        1.273392  345613.340048  0.020473   
7  pipe2_tj_350_tf_350     350.0        1.273392  345613.340048  0.020473   

         dp   dp_calc  
0  1.737575  1.737575  
1  1.737575  1.737575  
2  1.686985  1.686985  
3  1.659683  1.659683  
4  1.686985  1.686985  
5  1.737575  1.737575  
6  1.659683  1.659683  
7  1.659683  1.659683  

The output demonstrates:

• v_mean_m_per_s does not depend on the temperatures and is calculated with the density at 273.15 K.
• reynolds and in consequence lambda depend on the tfluid_k at the junctions and from the external grid.
  • reynolds for pipe1 has three different values: When both temperatures are low => lowest value; when one of the mentioned temperatures is high => medium value; when both temperature are high => highest value.
  • reynolds for pipe2 has only two values, depending on tfluid_k at the junctions and ignores the actual fluid temperature totally.
    I guess, that this is an effect of using tfluid_k and t_k at external grids for the calculation of the viscosity (instead of the actual temperature of the fluid in the pipes).
• The last column in the output dataframe 'dp_calc' shows, that these unexpected values for velocity and lambda are actually used for the pressure drop calculation.

Expected behavior
Use the actual temperatures of the fluid in the pipe to calculate the fluid properties and results, that depend on those (velocity, reynolds, lambda, pressure drop).

Python environment (please complete the following information):

• OS: Windows 10
• pandapipes version 0.10.0
• Note that in pandapipes version 0.9.0 the output of the above example is somewhat different: (probably related to Getting different velocities in versions 0.9.0 and 0.10.0 #625)
  Here velocity and reynolds only depend on tfluid_k at junctions: (which is somewhat more consistent than the pp 0.10.0 behavior)

All flow velocities equal: False
Calculated densities from v_mean: {996.49, 973.62}
Actual density at 273.15 K: 999.8801666666666
Actual density at 300 K: 996.49
Actual density at 350 K: 973.62
                  name  t_from_k  v_mean_m_per_s       reynolds    lambda  \
0  pipe1_tj_300_tf_300     300.0        1.277724  149231.076504  0.021413
1  pipe2_tj_300_tf_300     300.0        1.277724  149231.076504  0.021413
2  pipe1_tj_350_tf_300     300.0        1.307738  345613.340048  0.020463
3  pipe2_tj_350_tf_300     300.0        1.307738  345613.340048  0.020463
4  pipe1_tj_300_tf_350     350.0        1.277724  149231.076504  0.021413
5  pipe2_tj_300_tf_350     350.0        1.277724  149231.076504  0.021413
6  pipe1_tj_350_tf_350     350.0        1.307738  345613.340048  0.020463
7  pipe2_tj_350_tf_350     350.0        1.307738  345613.340048  0.020463

         dp   dp_calc
0  1.741814  1.741814
1  1.741814  1.741814
2  1.703636  1.703636
3  1.703636  1.703636
4  1.741814  1.741814
5  1.741814  1.741814
6  1.703636  1.703636
7  1.703636  1.703636

Additional context
This issue is related to #384, and illustrates, that the fully separated hydraulic and thermal calculation as implemented now leads to confusing results.

[dlohmeier]

On the develop branch, we implemented a new calculation mode "bidirectional" which would probably lead to more expected results. If you would like to test it with your example, we gladly use your feedback to make further improvements. I am not quite sure how to interpret the "real" pipe temperature, as temperature is hard to grasp for a component along which it changes. It will always be some intermediate value and must be derived from specific locations, which are usually the network nodes, in my opinion. But I would be interested in further ideas.

[dlohmeier]

I just found that there was a bug in the result extraction for incompressible media. I opened pull request #640 as a fix. It will be on develop soon, I will also clarify how we can release it rather soon.

[j-zipplies]

I totally agree, that "the" temperature in a pipe does not exist.
In district heating the temperature difference between inlet and outlet should typically be <1K. Thus, it will be a very good approximation to use the mean temperature between inlet and outlet of the pipe (not mixing temperature at the outlet node) for the calculation of fluid properties, as you propose. For even more accuracy - e.g. in long pipes and during low load periods, where temperature drops might be higher - the pipe could be segmented and fluid properties calculated segment-wise.
Another case are the components with larger temperature difference: heat exchangers or circ pumps. Here, density and velocity are not important, as these components don't calculate a pressure drop based on hydraulic flow calculations (as far as I understand them). However, the (temperature dependent) heat capacity is important in those components. High accuracy could be achieved by splitting the temperature drop into smaller portions and calculate cp as a mean value. E.g. a heat exchanger with 80 °C Inlet, 50 °C outlet => cp_mean = mean( cp(75°C), cp(65°C), cp(55°C) ). I don't know, how this is done currently.

[j-zipplies]

If I use your bugfix from #640 and the mode "bidirectional" that you mentioned, I actually get results, that look more like I expected:

• All properties and thus the results for fluid velocity and pressure drop depend on the actual temperature of the fluid, such that the first four pipes (300 K) and the second four pipes (350 K) have the same results.
• However, the check of calculating dp_calc = lambda * l / d * 0.5 * rho * v² yields different results (that I think are the correct ones, I compared to another pressure drop calculation tool https://www.schweizer-fn.de/berechnung/stroemung/druckverlust_rohr_rech.php). I noted, that the factor between dp / dp_calc is density(t_from_k) / density(273.15). So, I guess that at some point in the pressure drop calculation still the density at 273.15 K is used.

Calculated densities from v_mean: {996.49, 973.62}
Actual density at 273.15 K: 999.8801666666666
Actual density at 300 K: 996.49
Actual density at 350 K: 973.62
                  name  t_from_k  v_mean_m_per_s       reynolds    lambda  \
0  pipe1_tj_300_tf_300     300.0        1.277724  149231.076504  0.021434
1  pipe2_tj_300_tf_300     300.0        1.277724  149231.076504  0.021434
2  pipe1_tj_350_tf_300     300.0        1.277724  149231.076504  0.021434
3  pipe2_tj_350_tf_300     300.0        1.277724  149231.076504  0.021434
4  pipe1_tj_300_tf_350     350.0        1.307738  345613.340048  0.020473
5  pipe2_tj_300_tf_350     350.0        1.307738  345613.340048  0.020473
6  pipe1_tj_350_tf_350     350.0        1.307738  345613.340048  0.020473
7  pipe2_tj_350_tf_350     350.0        1.307738  345613.340048  0.020473

         dp   dp_calc
0  1.737575  1.743487
1  1.737575  1.743487
2  1.737575  1.743487
3  1.737575  1.743487
4  1.659683  1.704448
5  1.659683  1.704448
6  1.659683  1.704448
7  1.659683  1.704448

[dlohmeier]

@SimonRubenDrauz can you check the following line of code again:

pandapipes/src/pandapipes/pf/derivative_toolbox.py

Line 22 in 5b7e42c

const_term = np.divide(1, branch_pit[:, AREA] ** 2 * rho_n * P_CONVERSION * 2)

--> there we use density at norm state, and I have in mind that there was some purpose to it. Can we review this together? Thanks a lot @j-zipplies for the helpful feedback! I will try to add the required changes to the fix PR.

[dlohmeier]

With respect to the heat capacity, I could imagine that it should be part of the fluid property itself to find the mean value within a certain range, as the behavior can change a lot. So, if the dependency is e.g. quadratic, the property will have to find the integral for the given range and divide it by the range to retrieve the mean value. In case of interpolated points, the points within the range have to be evaluated and weighted by the respective distances.

[dlohmeier]

Hi @j-zipplies, I corrected the bug and added your snippet as a test so that this issue will not evolve again in future versions (a143789). This bug fix will be released soon. Thanks a lot for your support!