Add examples.

examples/script_design_capacity.py

+#******************************************************************************
+#******************************************************************************
+
+# TIP: this script requires the data found in the hyhetra-pipedata and hyhetra-
+# -fluiddata repositories
+
+pipedata_files = ['pipes/isoplus_singlepipes_s1.csv']
+        
+fluiddata_file = 'fluids/incropera2006_saturated_water.csv'
+
+# other files
+
+file_best2018 = 'papers/Best2018_design_flow_speeds.csv'
+
+file_sven2013 = 'papers/Sven2013_design_flow_speeds.csv'
+
+file_nussbaumer2016 = 'papers/Nussbaumer2016_design_flow_speeds.csv'
+
+file_euroheat2021 = 'papers/EuroHeat2021_design_flow_speeds.csv'
+
+# import libraries
+import numpy as np
+import matplotlib.pyplot as plt
+# topupheat
+import topupheat.pipes.fic as fic
+from topupheat.pipes.single import StandardisedPipe, StandardisedPipeDatabase
+import topupheat.pipes.utils as utils
+from topupheat.common import formulas as core
+
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+
+# specifications
+
+# maximum specific pressure drop
+
+list_max_specific_pressure_loss = [100, 200, 200, 200, 450] # Pa/m
+
+list_pipe_roughness = [0.1/1000, 0.1/1000, 0.05/1000, 0.01/1000, 0.1/1000]
+
+# district heating details
+
+dh_flow_temperature = 273.15+80 # K
+
+dh_return_temperature = 273.15+50 # K
+
+# temperature fluid bulk
+
+temperature_fluid_bulk = 0.5*(dh_return_temperature+dh_flow_temperature)
+
+# ground temperature
+
+ground_temperature = 273.15+10 # K
+
+# pipe absolute/effective roughness
+
+pipe_e_eff = 0.01/1000 # m
+
+#******************************************************************************
+
+# specify the pipes that should be considered using (DN,S) tuples
+
+list_pipe_tuples = [(20, 1, 'isoplus', 'DRE-20-STD'),
+                    (25, 1, 'isoplus', 'DRE-25-STD'),
+                    (32, 1, 'isoplus', 'DRE-32-STD'),
+                    (40, 1, 'isoplus', 'DRE-40-STD'),
+                    (50, 1, 'isoplus', 'DRE-50-STD'),
+                    (65, 1, 'isoplus', 'DRE-65-STD'),
+                    (80, 1, 'isoplus', 'DRE-80-STD'),
+                    (100, 1, 'isoplus', 'DRE-100-ST'),
+                    (125, 1, 'isoplus', 'DRE-125-ST'),
+                    (150, 1, 'isoplus', 'DRE-150-ST'),
+                    (200, 1, 'isoplus', 'DRE-200-ST'),
+                    (250, 1, 'isoplus', 'DRE-250-ST'),
+                    (300, 1, 'isoplus', 'DRE-300-ST'),
+                    (350, 1, 'isoplus', 'DRE-350-ST'),
+                    (400, 1, 'isoplus', 'DRE-400-ST'),
+                    (450, 1, 'isoplus', 'DRE-450-ST'),
+                    (500, 1, 'isoplus', 'DRE-500-ST'),
+                    (600, 1, 'isoplus', 'DRE-600-ST'),
+                    (700, 1, 'isoplus', 'DRE-700-ST'),
+                    (800, 1, 'isoplus', 'DRE-800-ST'),
+                    (900, 1, 'isoplus', 'DRE-900-ST'),
+                    (1000, 1, 'isoplus', 'DRE-1000-S')]
+
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+
+# pipe data
+
+pipedb = StandardisedPipeDatabase(source=pipedata_files)
+
+# water data
+
+waterdb = fic.FluidDatabase(fluid='water', phase='l', source=fluiddata_file)
+
+#******************************************************************************
+#******************************************************************************
+
+# create the objects for the fluids in the flow and return pipes
+
+dh_flow = fic.Fluid(
+    phase='l',
+    temperature=dh_flow_temperature,
+    pressure=1e5,
+    db=waterdb)
+
+dh_return = fic.Fluid(
+    phase='l',
+    temperature=dh_return_temperature,
+    pressure=1e5,
+    db=waterdb)
+
+fluid_bulk = fic.Fluid(
+    phase='l',
+    temperature=temperature_fluid_bulk,
+    pressure=1e5,
+    db=waterdb)
+            
+#******************************************************************************
+#******************************************************************************
+
+# define the pipes under consideration for task 2
+
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+
+# 4) compute the maximum heat flow rate for a given maximum specific pressure drop
+# as a function of pipe diameters
+                   
+#******************************************************************************
+
+# dictionaries
+
+dict_heat_transfer_rate = {}
+
+dict_friction_factor = {}
+
+dict_fluid_speed = {}
+
+dict_reynolds_number = {}
+
+dict_specific_pressure_loss = {}
+
+# for each specific pressure loss level
+
+for j, max_specific_pressure_loss in enumerate(list_max_specific_pressure_loss):
+
+    # declare the final and intermediate variables
+    
+    list_heat_transfer_rate = []
+    
+    list_friction_factor = []
+    
+    list_fluid_speed = []
+    
+    list_reynolds_number = []
+    
+    list_specific_pressure_loss = []
+    
+    temp_fluid_speed, temp_reynolds_number, temp_friction_factor = 0, 0, 0
+    
+    # pipe object list
+    
+    list_dh_pipes = [
+        StandardisedPipe(
+            pipe_tuple=pipe_tuple,
+            e_eff=list_pipe_roughness[j],
+            db=pipedb)
+        for pipe_tuple in list_pipe_tuples]
+        
+    # for each pipe under consideration
+    
+    for i, pipe in enumerate(list_dh_pipes):
+                                            
+        #**********************************************************************
+        
+        # get an initial solution
+        # >> should be valid and as close to the actual solution as possible
+        # >> one option is to use explicit methods
+        
+        if i == 0:
+            
+            x0_fluid_speed = None
+            x0_reynolds_number = None
+            x0_friction_factor = None
+            
+        else:
+            
+            # use previous solution
+            
+            (x0_fluid_speed,
+             x0_reynolds_number,
+             x0_friction_factor) = (list_fluid_speed[i-1],
+                                    list_reynolds_number[i-1],
+                                    list_friction_factor[i-1])
+        
+        # flow regime
+        
+        # compute the flow speed that leads to the specified pressure loss
+        
+        temp_fluid_speed, temp_reynolds_number, temp_friction_factor = (
+            utils.find_specific_pressure_loss_matching_flow(
+                pipe, 
+                fluid_bulk,
+                max_specific_pressure_loss,
+                dff_model=fic.DFF_COLEBROOK_WHITE,
+                #dff_model=fic.DFF_HAALAND,
+                #dff_model=fic.DFF_NIKURADSE,
+                #dff_model=fic.DFF_PETHUKOV,
+                #dff_model=fic.DFF_BLASIUS,
+                #dff_model=fic.DFF_PRANDTL_NIKURADSE_KARMAN,
+                x0_friction_factor=x0_friction_factor,
+                x0_reynolds_number=x0_reynolds_number,
+                x0_fluid_speed=x0_fluid_speed)
+            )  
+        
+        temp_specific_pressure_loss = fic.SpecificPressureLossInPipe(
+            fluid_bulk.mass_density, 
+            pipe.d_int, 
+            temp_fluid_speed, 
+            temp_friction_factor)
+        
+        # store the results
+        
+        list_fluid_speed.append(temp_fluid_speed)
+        
+        list_reynolds_number.append(temp_reynolds_number)
+        
+        list_friction_factor.append(temp_friction_factor)
+        
+        list_specific_pressure_loss.append(temp_specific_pressure_loss)
+        
+        # compute the heat flow rate
+        
+        temp_heat_flow_rate = fic.HeatTransferRateInPipe(
+            temp_fluid_speed, 
+            dh_flow, 
+            dh_return, 
+            pipe,
+            fluid_bulk=fluid_bulk
+            )
+        
+        # store the result
+        
+        list_heat_transfer_rate.append(temp_heat_flow_rate)
+        
+    #**************************************************************************
+        
+    # add lists to the dicts
+    
+    dict_specific_pressure_loss.update(
+        {(max_specific_pressure_loss,list_pipe_roughness[j]):
+         list_specific_pressure_loss})
+        
+    dict_fluid_speed.update(
+        {(max_specific_pressure_loss,list_pipe_roughness[j]):
+         list_fluid_speed})
+    
+    dict_reynolds_number.update(
+        {(max_specific_pressure_loss,list_pipe_roughness[j]):
+         list_reynolds_number})
+    
+    dict_friction_factor.update(
+        {(max_specific_pressure_loss,list_pipe_roughness[j]):
+         list_friction_factor})
+    
+    dict_heat_transfer_rate.update(
+        {(max_specific_pressure_loss,list_pipe_roughness[j]):
+         list_heat_transfer_rate})
+    
+    #**************************************************************************
+    
+    # print('******************************************************************')
+    # print('******************************************************************')
+    # print('******************************************************************')
+    # print('******************************************************************')
+        
+#******************************************************************************
+
+isoplus_diameters = [20,
+                     25,
+                     32,
+                     40,
+                     50,
+                     65,
+                     80,
+                     100,
+                     125,
+                     150,
+                     200,
+                     250,
+                     300,
+                     350,
+                     400,
+                     450,
+                     500,
+                     600,
+                     700,
+                     800,
+                     900,
+                     1000]
+
+isoplus_mass_flow_rate_min_t_hour = [0.4,
+                                     0.8,
+                                     1.7,
+                                     2.5,
+                                     4.7,
+                                     9.3,
+                                     14.5,
+                                     28.5,
+                                     50.0,
+                                     82.0,
+                                     167.0,
+                                     300,
+                                     472,
+                                     610,
+                                     862,
+                                     1180,
+                                     1570,
+                                     2520,
+                                     3770,
+                                     5390,
+                                     7400,
+                                     9200]
+
+isoplus_mass_flow_rate_max_t_hour = [0.5,
+                                     1.0,
+                                     2.0,
+                                     3.0,
+                                     5.5,
+                                     11.0,
+                                     16.5,
+                                     33.0,
+                                     58.0,
+                                     95.0,
+                                     193.0,
+                                     348,
+                                     547,
+                                     705,
+                                     1000,
+                                     1370,
+                                     1820,
+                                     2920,
+                                     4370,
+                                     6240,
+                                     9500,
+                                     np.nan]
+
+isoplus_fluid_speed_min = [
+    core.MeanFluidSpeedViaVolumetricFlowRate(
+        core.VolumetricFlowRateFromMassFlowRate(
+            m*0.27777778, #t/h to kg/s
+            fluid_bulk.mass_density
+            ), 
+        fic.Pipe.CrossSectionalArea(fic.Pipe, isoplus_diameters[i]/1e3)
+        )
+    for i, m in enumerate(isoplus_mass_flow_rate_min_t_hour)]
+
+isoplus_fluid_speed_max = [
+    core.MeanFluidSpeedViaVolumetricFlowRate(
+        core.VolumetricFlowRateFromMassFlowRate(
+            m*0.27777778, #t/h to kg/s
+            fluid_bulk.mass_density
+            ), 
+        fic.Pipe.CrossSectionalArea(fic.Pipe, isoplus_diameters[i]/1e3)
+        )
+    for i, m in enumerate(isoplus_mass_flow_rate_max_t_hour)]
+
+# get data from Best et al. (2018)
+
+data_best2018 = np.genfromtxt(file_best2018,
+                              #dtype=(float,float,float,float),
+                              names=True,
+                              # skip_header=1,
+                              # skip_footer=0,
+                              delimiter=',')
+
+data_sven2013 = np.genfromtxt(file_sven2013,
+                              #dtype=(float,float,float,float),
+                              names=True,
+                              # skip_header=1,
+                              # skip_footer=0,
+                              delimiter=',')
+
+data_nussbaumer2016 = np.genfromtxt(
+    file_nussbaumer2016,
+    #dtype=(float,float,float,float),
+    names=True,
+    # skip_header=1,
+    # skip_footer=0,
+    delimiter=',')
+
+data_euroheat2021 = np.genfromtxt(
+    file_euroheat2021,
+    #dtype=(float,float,float,float),
+    names=True,
+    # skip_header=1,
+    # skip_footer=0,
+    delimiter=',')
+
+# list_pipe_diameters = [t[0] for t in list_pipe_tuples]
+
+list_pipe_diameters = [pipe.d_int*1000 for pipe in list_dh_pipes]
+
+# plot the max heat flow rate as function of pipe diameter (1st pos. in tuple)
+
+fig, ax = plt.subplots()
+
+for i, max_specific_pressure_loss in enumerate(list_max_specific_pressure_loss):
+
+    ax.plot(list_pipe_diameters,
+            [-y/1e6 for y in dict_heat_transfer_rate[
+                (max_specific_pressure_loss,list_pipe_roughness[i])]],
+            label='MSP='+str(max_specific_pressure_loss)+' Pa/m')
+
+ax.set(xlabel='Pipe internal diameter [mm]', 
+       ylabel='Heat Transfer Rate [MW]',
+       title='Maximum heat capacity')
+
+ax.legend()
+
+ax.grid()
+
+# fig.savefig("test.png")
+plt.show()
+    
+#******************************************************************************    
+
+# plot the max fluid speed as a function of pipe diameter (1st pos. in tuple)
+
+fig, ax = plt.subplots()
+
+fig.set_size_inches(12,8) # (15,10) (12,8) # (3,2) (3*2,2*2) (3*3,2*3)
+
+list_linestyle_our = ['rx-','bx-','gx-','cx-','yx-']
+
+list_linestyle_euroheat2021 = ['rs','bs','gs','cs','ys']
+
+list_linestyle_best2018 = ['r>--','b>--','g>--']
+
+list_linestyle_nussbaumer2016 = ['r<-.','b<-.','g<-.']
+
+published_results_labels_euroheat2021 = ['100_Pa_m_01mm',
+                                         '200_Pa_m_01mm',
+                                         '200_Pa_m_005mm',
+                                         '200_Pa_m_001mm',
+                                         '450_Pa_m_01mm']
+
+published_results_labels_best2018 = ['100_Pa_m','200_Pa_m','300_Pa_m']
+
+published_results_labels_nussbaumer2016 = ['100_Pa_m','200_Pa_m','300_Pa_m']
+
+legend_euroheat2021 = ['100 Pa/m, e=0.1mm',
+                       '200 Pa/m, e=0.1mm',
+                       '200 Pa/m, e=0.05mm',
+                       '200 Pa/m, e=0.01mm',
+                       '450 Pa/m, e=0.1mm']
+
+legend_best2018 = ['100 Pa/m, e=0.01mm',
+                   '200 Pa/m, e=0.01mm',
+                   '300 Pa/m, e=0.01mm']
+
+legend_nussbaumer2016 = ['100 Pa/m, e=0.01mm',
+                         '200 Pa/m, e=0.01mm',
+                         '300 Pa/m, e=0.01mm']
+
+# x,100_Pa_m_01mm,200_Pa_m_01mm,200_Pa_m_005mm,200_Pa_m_001mm,450_Pa_m_01mm
+
+for i, max_specific_pressure_loss in enumerate(list_max_specific_pressure_loss):
+
+    ax.plot(list_pipe_diameters,
+            [fluid_speed for fluid_speed in dict_fluid_speed[
+                (max_specific_pressure_loss,list_pipe_roughness[i])]],
+            #marker='s', linestyle='-',
+            list_linestyle_our[i],
+            label='topupheat, '+str(max_specific_pressure_loss)+' Pa/m, e='+str(list_pipe_roughness[i]*1e3)+' mm')
+
+for i, label in enumerate(published_results_labels_euroheat2021):
+            
+    ax.plot(data_euroheat2021['x']*1e3,
+            data_euroheat2021[label],
+            #marker='x', linestyle='-.',
+            list_linestyle_euroheat2021[i],
+            label='EuroHeat&Power (2021), '+str(legend_euroheat2021[i]))
+
+# for i, label in enumerate(published_results_labels_best2018):
+            
+#     ax.plot(data_best2018['x']*1e3,
+#             data_best2018[label],
+#             #marker='x', linestyle='-.',
+#             list_linestyle_best2018[i],
+#             label='Best et al. (2018), '+str(legend_best2018[i]))
+
+# for i, label in enumerate(published_results_labels_nussbaumer2016):
+            
+#     ax.plot(data_nussbaumer2016['x']*1e3,
+#             data_nussbaumer2016[label],
+#             #marker='x', linestyle='-.',
+#             list_linestyle_nussbaumer2016[i],
+#             label='Nussbaumer et al. (2016), '+str(legend_nussbaumer2016[i]))
+
+# ax.plot(data_sven2013['x']*1e3,
+#         data_sven2013['100_Pa_m'],
+#         #marker='x', linestyle='-.',
+#         'kd-.',
+#         label='Svend and Sven (2013), 100 Pa/m')
+
+#******************************************************************************
+ax.set(xlim=(0, 1000),ylim=(0, 9))
+ax.set(xlabel='Pipe internal diameter [mm]', 
+       ylabel='Fluid speed [m/s]',
+       title='Maximum mean fluid speed for a given specific pressure loss')
+
+ax.legend()
+
+ax.grid()
+
+# fig.savefig("test.png")
+plt.show()
+    
+#******************************************************************************    
+
+# plot the specific pressure loss as a function of pipe diameter
+
+fig, ax = plt.subplots()
+
+for i, max_specific_pressure_loss in enumerate(list_max_specific_pressure_loss):
+
+    ax.plot(list_pipe_diameters,
+            [spl for spl in dict_specific_pressure_loss[
+                (max_specific_pressure_loss,list_pipe_roughness[i])]],
+            label='MSP='+str(max_specific_pressure_loss)+' Pa/m')
+
+    # ax.plot(data_best2018['x']*1e3,
+    #         data_best2018[published_results_labels[i]],
+    #         #marker='o', linestyle='--',
+    #         list_linestyle_best2018[i],
+    #         label='Best et al. (2018), '+str(max_specific_pressure_loss)+' Pa/m')
+    
+    # ax.plot(data_nussbaumer2016['x']*1e3,
+    #         data_nussbaumer2016[published_results_labels[i]],
+    #         #marker='x', linestyle='-.',
+    #         list_linestyle_nussbaumer2016[i],
+    #         label='Nussbaumer et al. (2016), '+str(max_specific_pressure_loss)+' Pa/m')
+
+    
+ax.set(xlabel='Pipe internal diameter [mm]', 
+       ylabel='Specific pressure loss [Pa/m]',
+       title='...')
+
+ax.legend()
+
+ax.grid()
+
+# fig.savefig("test.png")
+plt.show()
+
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
\ No newline at end of file

examples/script_moody_diagram.py

+#******************************************************************************
+#******************************************************************************
+
+# TIP: this script requires the data found in the hyhetra-pipedata and hyhetra-
+# -fluiddata repositories
+
+pipedata_files = ['pipes/isoplus_singlepipes_s1.csv']
+
+fluiddata_file = 'fluids/incropera2006_saturated_water.csv'
+
+# import libraries
+import numpy as np
+import matplotlib.pyplot as plt
+import topupheat.pipes.fic as fic
+from topupheat.pipes.single import StandardisedPipe, StandardisedPipeDatabase
+
+#******************************************************************************
+#******************************************************************************
+
+# specify the pipes that should be considered 
+
+list_pipe_tuples = [(20, 1, 'isoplus', 'DRE-20-STD'),
+                    (25, 1, 'isoplus', 'DRE-25-STD'),
+                    (32, 1, 'isoplus', 'DRE-32-STD'),
+                    (40, 1, 'isoplus', 'DRE-40-STD'),
+                    (50, 1, 'isoplus', 'DRE-50-STD'),
+                    (65, 1, 'isoplus', 'DRE-65-STD'),
+                    (80, 1, 'isoplus', 'DRE-80-STD'),
+                    (100, 1, 'isoplus', 'DRE-100-ST'),
+                    (125, 1, 'isoplus', 'DRE-125-ST'),
+                    (150, 1, 'isoplus', 'DRE-150-ST'),
+                    (200, 1, 'isoplus', 'DRE-200-ST')]
+
+# relative pipe roughness
+
+list_relative_pipe_roughness = [
+    5e-2,
+    2.5e-2,
+    1e-2,
+    5e-3,
+    2e-3,
+    1e-3,
+    5e-4,
+    2e-4,
+    1e-4,
+    5e-5,
+    2e-5,
+    1e-5,
+    0
+    ]
+
+labels_relative_pipe_roughness = [
+    '5e-2',
+    '2.5e-2',
+    '1e-2',
+    '5e-3',
+    '2e-3',
+    '1e-3',
+    '5e-4',
+    '2e-4',
+    '1e-4',
+    '5e-5',
+    '2e-5',
+    '1e-5',
+    '0']
+
+# specify the pipe that is to be considered individually
+
+selected_pipe_index = 5
+
+# moody diagram
+
+number_points_laminar = 5
+
+number_points_turbulent = 100
+
+#******************************************************************************
+#******************************************************************************
+
+# pipe data
+
+pipedb = StandardisedPipeDatabase(source=pipedata_files)
+
+# water data
+
+waterdb = fic.FluidDatabase(fluid='water', phase='l', source=fluiddata_file)
+
+#******************************************************************************
+#******************************************************************************
+
+# 1) plot the moody diagram
+
+#******************************************************************************
+
+# define the reynolds numbers' vector
+
+list_reynolds_number = (
+    np.logspace(0,
+                8, 
+                number_points_turbulent+number_points_laminar)
+    )
+
+# define the models to be used
+
+list_dff_models = [fic.DFF_COLEBROOK_WHITE]
+
+#******************************************************************************
+#******************************************************************************
+
+# variables
+
+dict_friction_factor = {}
+
+dict_pipe = {}
+
+# for each roughness value
+
+for relative_pipe_roughness in list_relative_pipe_roughness:
+    
+    # create the pipe instance
+    
+    pipe = StandardisedPipe(
+        pipe_tuple=list_pipe_tuples[selected_pipe_index],
+        e_rel=relative_pipe_roughness,
+        db=pipedb)
+    
+    dict_pipe[relative_pipe_roughness] = pipe
+
+    # for each dff model
+
+    for model in list_dff_models:
+    
+        # get the friction factors
+        
+        list_friction_factor = [fic.DarcyFrictionFactor(
+            re, 
+            pipe,
+            model=model) 
+            for re in list_reynolds_number]
+        
+        # store it
+        
+        dict_friction_factor[
+            (model, relative_pipe_roughness)] = list_friction_factor
+            
+        #**********************************************************************
+        #**********************************************************************
+        
+    #**************************************************************************
+    #**************************************************************************
+
+#******************************************************************************
+#******************************************************************************
+
+list_linestyle_our = ['ro-','go-','bx-',
+                      'rv-','g^-','bd-',
+                      'r<-','g>-','b.-',
+                      'ro--','go--','bx--',
+                      'rv--','g^--','bd--',
+                      'r<--','g>--','b.--']
+
+# plot
+
+fig, ax = plt.subplots()
+
+fig.set_size_inches(12,8)
+    
+for j, model in enumerate(list_dff_models):
+
+    for i, relative_pipe_roughness in enumerate(list_relative_pipe_roughness):
+
+        # plot them using a log scale
+        
+        linestyle_index = i+(j)*len(list_relative_pipe_roughness)
+        
+        ax.loglog(list_reynolds_number, 
+                  dict_friction_factor[(model, relative_pipe_roughness)],
+                  list_linestyle_our[linestyle_index],
+                  #label=model+', e/D='+labels_relative_pipe_roughness[i],
+                  label='CWE, e/D='+labels_relative_pipe_roughness[i])
+        
+        #**********************************************************************
+        
+    #**************************************************************************
+        
+#******************************************************************************
+#******************************************************************************
+    
+ax.set(xlabel='Reynolds number [-]', 
+       ylabel='Friction factor [-]',
+       title='Moody diagram')
+
+# plt.autoscale(enable=True, axis='x', tight=True)
+
+ax.grid(which='both')
+
+ax.legend() 
+
+# left, right = xlim()  # return the current xlim
+
+plt.xlim((0.6e3, 1e8))   # set the xlim to left, right
+plt.ylim((0.006, 0.1))   # set the xlim to left, right
+
+# fig.savefig("test.png")
+plt.show()
+    
+#******************************************************************************
+#******************************************************************************
\ No newline at end of file

examples/script_nikuradse_experiments.py

+#******************************************************************************
+#******************************************************************************
+
+# TIP: this script requires the data found in the hyhetra-pipedata and hyhetra-
+# -fluiddata repositories
+
+pipedata_files = ['pipes/isoplus_singlepipes_s1.csv']
+
+fluiddata_file = 'fluids/incropera2006_saturated_water.csv'
+
+# Note: this script uses data published by Li and Huai (2016) about Nikuradse's
+# experiments. DOI:  https://doi.org/10.1371/journal.pone.0154408
+
+# The xls file included in the paper as an appendix should then be converted to
+# csv format and its two first lines (headers) replaced by the following one:
+# A_run,A_f,A_1/e,A_Re,B_run,B_f,B_1/e,B_Re,C_run,C_f,C_1/e,C_Re,D_run,D_f,D_1/e,D_Re,E_run,E_f,E_1/e,E_Re,F_run,F_f,F_1/e,F_Re    
+# How? copy the line above (minus the leading '# ') and paste it using a text editor
+
+nikuradse_file = 'papers/Li2016_edited.csv'
+
+# import libraries
+import numpy as np
+import matplotlib.pyplot as plt
+import topupheat.pipes.fic as fic
+from topupheat.pipes.single import StandardisedPipe, StandardisedPipeDatabase
+
+#******************************************************************************
+#******************************************************************************
+
+# specify the pipes that should be considered using (DN,S) tuples
+
+# 20,25,32,40,50,65,80,100,125,150,200
+
+list_pipe_tuples = [(20, 1, 'isoplus', 'DRE-20-STD'),
+                    (25, 1, 'isoplus', 'DRE-25-STD'),
+                    (32, 1, 'isoplus', 'DRE-32-STD'),
+                    (40, 1, 'isoplus', 'DRE-40-STD'),
+                    (50, 1, 'isoplus', 'DRE-50-STD'),
+                    (65, 1, 'isoplus', 'DRE-65-STD'),
+                    (80, 1, 'isoplus', 'DRE-80-STD'),
+                    (100, 1, 'isoplus', 'DRE-100-ST'),
+                    (125, 1, 'isoplus', 'DRE-125-ST'),
+                    (150, 1, 'isoplus', 'DRE-150-ST'),
+                    (200, 1, 'isoplus', 'DRE-200-ST')]
+
+# relative pipe roughness
+
+list_relative_pipe_roughness = [
+    0.5/507,
+    0.5/252,
+    0.5/126,
+    0.5/60,
+    0.5/30.6,
+    0.5/15]
+
+labels_relative_pipe_roughness = [
+    '1/1014',
+    '1/504',
+    '1/252',
+    '1/120',
+    '1/61.2',
+    '1/30']
+
+# specify the pipe that is to be considered individually
+
+selected_pipe_index = 5
+
+# moody diagram
+
+number_points_laminar = 5
+
+number_points_turbulent = 100
+
+#******************************************************************************
+#******************************************************************************
+
+# pipe data
+
+pipedb = StandardisedPipeDatabase(source=pipedata_files)
+
+# water data
+
+waterdb = fic.FluidDatabase(fluid='water', phase='l', source=fluiddata_file)
+
+#******************************************************************************
+#******************************************************************************
+
+nikuradse_data = np.genfromtxt(nikuradse_file,
+                               #dtype=(float,float,float,float),
+                               names=True,
+                               #skip_header=0,
+                               # skip_header=1,
+                               # skip_footer=0,
+                               delimiter=',')
+
+nikuradse_experiments = ['A','B','C','D','E','F']
+
+nikuradse_data_length = [48, 38, 71, 72, 63, 70]
+
+nikuradse_data_f = {nikuradse_experiments[i]: 
+                    nikuradse_data[value+'_f'][0:nikuradse_data_length[i]]
+                    for i, value in enumerate(nikuradse_experiments)}
+
+nikuradse_data_re = {nikuradse_experiments[i]: 
+                     nikuradse_data[value+'_Re'][0:nikuradse_data_length[i]]
+                     for i, value in enumerate(nikuradse_experiments)}
+
+nikuradse_data_1e = {nikuradse_experiments[i]: 
+                     nikuradse_data[value+'_1e'][1]
+                     for i, value in enumerate(nikuradse_experiments)}
+
+#******************************************************************************
+#******************************************************************************
+
+# 1) plot the moody diagram
+
+#******************************************************************************
+
+# define the reynolds numbers' vector
+
+list_reynolds_number = (
+    np.logspace(0,
+                8, 
+                number_points_turbulent+number_points_laminar)
+    )
+
+# define the models to be used
+
+list_dff_models = [fic.DFF_COLEBROOK_WHITE]
+
+#******************************************************************************
+#******************************************************************************
+
+# variables
+
+dict_friction_factor = {}
+
+dict_pipe = {}
+
+# for each roughness value
+
+for relative_pipe_roughness in list_relative_pipe_roughness:
+    
+    # create the pipe instance
+    
+    pipe = StandardisedPipe(
+        pipe_tuple=list_pipe_tuples[selected_pipe_index],
+        e_rel=relative_pipe_roughness,
+        db=pipedb)
+    
+    dict_pipe[relative_pipe_roughness] = pipe
+
+    # for each dff model
+
+    for model in list_dff_models:
+    
+        # get the friction factors
+        
+        list_friction_factor = [fic.DarcyFrictionFactor(
+            re, 
+            pipe,
+            model=model) 
+            for re in list_reynolds_number]
+        
+        # store it
+        
+        dict_friction_factor[
+            (model, relative_pipe_roughness)] = list_friction_factor
+            
+        #**********************************************************************
+        #**********************************************************************
+        
+    #**************************************************************************
+    #**************************************************************************
+
+#******************************************************************************
+#******************************************************************************
+
+list_linestyle_our = ['rs-','rd-','rx-',
+                      'r^-','r>-','r<-']
+
+list_linestyle_nikuradse = ['bs','bd','bx',
+                            'b^','b>','b<']
+
+# plot
+
+fig, ax = plt.subplots()
+
+fig.set_size_inches(12,8)
+    
+for j, model in enumerate(list_dff_models):
+
+    for i, relative_pipe_roughness in enumerate(list_relative_pipe_roughness):
+
+        # plot them using a log scale
+        
+        linestyle_index = i+(j)*len(list_relative_pipe_roughness)
+        
+        ax.loglog(list_reynolds_number, 
+                  dict_friction_factor[(model, relative_pipe_roughness)],
+                  list_linestyle_our[linestyle_index],
+                  #label=model+', e/D='+labels_relative_pipe_roughness[i],
+                  label='CWE, e/D='+labels_relative_pipe_roughness[i])
+        
+        #**********************************************************************
+        
+    #**************************************************************************
+    
+#******************************************************************************
+
+for e, exp in enumerate(nikuradse_experiments):
+    
+    ax.loglog(nikuradse_data_re[exp], 
+              nikuradse_data_f[exp],
+              list_linestyle_nikuradse[e],
+              label='N33, e/D=1/'+str(round(2*nikuradse_data_1e[exp])))
+    
+#******************************************************************************
+    
+ax.set(xlabel='Reynolds number [-]', 
+       ylabel='Friction factor [-]',
+       title='Colebrook-White equation (CWE) vs results from Nikuradse (1933)')
+
+# plt.autoscale(enable=True, axis='x', tight=True)
+
+ax.grid(which='both')
+
+ax.legend() 
+
+# left, right = xlim()  # return the current xlim
+
+# plt.xlim((1e3, 1e8))   # set the xlim to left, right
+plt.xlim((1e3, 1e6))   # set the xlim to left, right
+plt.ylim((1.5e-2, 8e-2))   # set the xlim to left, right
+
+# fig.savefig("test.png")
+plt.show()
+    
+#******************************************************************************
+#******************************************************************************
\ No newline at end of file