Preparation of cases

This notebooks aims to design the input files (json) for the cases presented in this paper. The idea is to automatically construct the input files fot the simulation, but considering the best initial estimations in order to guarantee the convergence.

With this the scenarios are two: short-term and long-term. The models are two: simplified and complete. The instances are three: preparation, simulation and simulationwithoutbiofilm. For more details regarding these definitions, please access the paper text.



In [1]:

    
# reading recquired modules
import json
from math import pi, ceil
from iapws.iapws97 import _Region4, _Region2, _Region1
from iapws._iapws import _ThCond, _Viscosity
from constructors import *
from scipy.optimize import root
import numpy as np
import os.path

# Constants
g = 9.81 # m/s2 - gravitational acceleration
pi = pi # pi

Cases



In [2]:

    
# Cases that where studied during the paper construction
v_cases = (0.9, 1.0, 1.15, 1.30) # approximate mean velocity for the heat exchanger tubes
kvap_cases = (7.077073, 7.077073, 7.077073, 7.077073) # steam mass flowrate (requires adjustment)
nodeA_P = [416246/4, 119929.5, 426067/4, 433048/4] # Estimated inlet network pressure
nodeA_w = [202.03, 224.43, 258.15, 291.82] # Estimate total cooling water flowrate
Tbfs = [312.5,310,309.7,308.0] # Estimate mean temperature of biofilm in heat exchanger
vbfs = [0.9, 1.0, 1.15, 1.30] # Estimate mean velocity of biofilm in heat exchanger


# Cases that will be considered in publication
required_cases = (1.0,)

Heat Exchanger

For more details regarding the geometry of the pipe network and the equipments, please access the paper text.

We are assuming a shell with a square cross section, single pass and stainless steel tubes (similar geometries can be found in http://www.doosanskodapower.com/), as presented in the following Figure (all rights to Doosan Škoda Power):



In [3]:

    
# Parameters
Do = 3/4*0.0254 # m - outter diameter 
Di = Do - 2*0.00165 # m - inner diameter 
L = 15 * 0.3048 # m - tube length
ep = 45e-6 # m - roughness 
kwall = 50.0 # steel thermal conductivity
external_volume = 20 # m3 - shell (steam) volume
KtubeInlet=0.5*0.9 # dischange coefficient for tube inlet
KtubeOutlet=0.5*0.9 # dischange coefficient for tube outlet
KnozzleInlet=1.1 # dischange coefficient for inlet nozzle
KnozzleOutlet=0.7 # dischange coefficient for outlet nozzle
At = 250 # m2 - Total heat exchange area of reference project (http://www.doosanskodapower.com/download/pdf/condensers.pdf (Plzeň project))

# Area per tube
Ahot = pi * Do * L # m2 - Tube heat exchange area
Acold = pi * Di * L # m2 - Tube heat exchange area
Ai = 0.25 * pi * Di**2 # m2 - Free flow area

# Tube bundle gemoetry
tubes_per_width = ceil((At/Acold)**0.5)
n_tubes = tubes_per_width ** 2
pitch = 2.5*Do # m inline
width = ((n_tubes ** 0.5) + 2) * pitch
tube_bundle = {
    "simplified": {
        "Ntubes_total": tubes_per_width ** 2,
        "Ntubes": 1,
        "Npipes": tubes_per_width ** 2,
    },
    "complete": {
        "Ntubes_total": n_tubes,
        "Ntubes": tubes_per_width,
        "Npipes": tubes_per_width,
    }
}

# Estimates (required for the preparation of the initial estimates)
# Surface temperature estimations
Ti0 = 312 # inner surface temperature at inlet
Ti1 = 312 # inner surface temperature at outlet
Ti = 0.5 * (Ti0+Ti1)
To0 = 314 # outter surface temperature at inlet
To1 = 314 # outter surface temperature at outlet
To = 0.5 * (To0+To1)

# Results
print("The tube inside diameter is {number:.{digits}f} m.".format(number=Di, digits=5))
print("The tube outside diameter is {number:.{digits}f} m.".format(number=Do, digits=5))
print("The tube length is {number:.{digits}f} m.".format(number=L, digits=2))
print("The tube roughness is {number:.{digits}f} m.".format(number=ep, digits=6))
print("The total heat exchange area (cold side) is {number:.{digits}f} m2.".format(number=At, digits=6))
print("The total number of tubes is {number:.{digits}f}.".format(number=n_tubes, digits=0))
print("The maximum number of inlined tubes is {number:.{digits}f}.".format(number=tubes_per_width, digits=0))
print("The shell width is {number:.{digits}f} m.".format(number=width, digits=2))









    



The tube inside diameter is 0.01575 m.
The tube outside diameter is 0.01905 m.
The tube length is 4.57 m.
The tube roughness is 0.000045 m.
The total heat exchange area (cold side) is 250.000000 m2.
The total number of tubes is 1156.
The maximum number of inlined tubes is 34.
The shell width is 1.71 m.

Exhaust Steam Quality

The exhaust steam quality of steam turbines may vary according to the designer. It is a common sense that the exhaust steam pressure shall be lower than atmospheric condition. According to Westinghouse steam turbine division, the turbine exhaust temperature should not exceed 175°F (70°C) for continuous operation or 250°F (121°C) for periods of about 15 minutes. So we are assuming the followinf Exhuast steam quality:



In [4]:

    
# Parameters
exhaust_steam_pressure = 10000 # Pa

# Estimated variable
estimated_steam_mass_flowrate = 4.392  # kg/s
estimated_steam_mass_flowrate = 7.078231821428575   # kg/s

# Steam properties
exhaust_steam_temperature = _Region4(exhaust_steam_pressure*1e-6, 0.)['T']
exhaust_steam_properties = _Region2(exhaust_steam_temperature,1e-6*exhaust_steam_pressure)
exhaust_steam_properties_liq = _Region1(exhaust_steam_temperature,1e-6*exhaust_steam_pressure)
exhaust_steam_density = 1/exhaust_steam_properties['v']
exhaust_steam_viscosity = _Viscosity(exhaust_steam_density, exhaust_steam_temperature)
exhaust_vaporization_heat = 1e3*exhaust_steam_properties['h'] - 1e3*exhaust_steam_properties_liq['h']

# Results
print("The exhaust steam pressure is equal to {0} Pa."\
      .format(exhaust_steam_pressure))
print("The exhaust steam temperature is equal to {number:.{digits}f}°C."\
      .format(number=exhaust_steam_temperature-273.15, digits=2))
print("The exhaust steam density is equal to {number:.{digits}f} kg/m3."\
      .format(number=exhaust_steam_density, digits=5))
print("The exhaust steam vaporization heat is equal to {number:.{digits}f} J/kg."\
      .format(number=exhaust_vaporization_heat, digits=0))









    



The exhaust steam pressure is equal to 10000 Pa.
The exhaust steam temperature is equal to 45.81°C.
The exhaust steam density is equal to 0.06816 kg/m3.
The exhaust steam vaporization heat is equal to 2392075 J/kg.

Pipe header



In [5]:

    
# Parameters
header_L = 50.
D_header  = 14 * 0.0254

Cooling water



In [6]:

    
# Parameters
water_inlet_temperature = 13.8 + 273.15 #K - cooling water inlet temperature (for the short-term scenario)
water_outlet_pressure = 100000 # Pa - cooling water outlet temperature

# Calculating properties
Tm = 300 # K - mean cooling water temperature
Pm = 100000 # Pa - mean cooling water pressure
water_properties_temperature = Tm
water_properties_pressure  = Pm
water_properties = _Region1(water_properties_temperature,1e-6*water_properties_pressure)
water_density = 1/water_properties['v']
water_entalphy = 1e3*water_properties['h']
water_heat_capacity = 1e3*water_properties['cp']
water_conductivity = _ThCond(water_density, water_properties_temperature)
water_viscosity = _Viscosity(water_density, water_properties_temperature)
water_prandtl = water_heat_capacity * water_viscosity / water_conductivity

# Results
print("Estimated mean cooling water temperature is {number:.{digits}f} K".format(number=Tm,digits=2))
print("The water density is {number:.{digits}f} kg/m3.".format(number=water_density, digits=2))
print("The water entalphy is {number:.{digits}f} J/kg.".format(number=water_entalphy, digits=2))
print("The water heat capacity is {number:.{digits}f} J/(kg*K).".format(number=water_heat_capacity, digits=2))
print("The water conductivity is {number:.{digits}f} W/(m*K).".format(number=water_conductivity, digits=3))
print("The water viscosity is {number:.{digits}f} (Pa*s).".format(number=water_viscosity, digits=5))
print("The water Prandtl is {number:.{digits}f}.".format(number=water_prandtl, digits=5))









    



Estimated mean cooling water temperature is 300.00 K
The water density is 996.56 kg/m3.
The water entalphy is 112663.82 J/kg.
The water heat capacity is 4181.10 J/(kg*K).
The water conductivity is 0.610 W/(m*K).
The water viscosity is 0.00085 (Pa*s).
The water Prandtl is 5.85657.

Biofilm



In [7]:

    
# Parameters
lagt = 1.0 # days - Lag time for biofilm formation
lagt_infinite = 1e10 # days - Forced Lag time to model the clean condenser
rhomf = 980.0 # kg/m3 - biofilm density for the simplified model
Ccell = 6e+7 # number of cells (NOT USED IN THIS MODEL)
Csubstrate = (1/10) * 20e-6 # mg/L - substrate concentration

Numerical Methods



In [8]:

    
# Parameters
Nelements = 10 # number of mesh points
simulation_parameters_default = {
"shortterm_preparation": {
    "reporting_interval": 10.,
    "time_horizon": 10.,
    "relative_tolerance":  1e-6,
    "MaxStep": 1,
    "MaxNumSteps": 1000000000},
"longterm_preparation": {
    "reporting_interval": 10.,
    "time_horizon": 10.,
    "relative_tolerance":  1e-6,
    "MaxStep": 1,
    "MaxNumSteps": 1000000000},
"shortterm_simulation": {
    "reporting_interval": 0.1,
    "time_horizon": 15,
    "relative_tolerance":  1e-6,
    "MaxStep": 1,
    "MaxNumSteps": 1000000000},
"longterm_simulation": {
    "reporting_interval": 1.0,
    "time_horizon": 365*2,
    "relative_tolerance":  1e-2,
    "MaxStep": 1,
    "MaxNumSteps": 1000000000},
}

Simplified model

Hydraulics Equations



In [9]:

    
def calculate_darcy(fD0, ep, D, Re):
    """Calculate the darcy friction factor"""
    def idarcy(ifD, ep, D, Re):
        return ifD + 2. * np.log10(ep / 3.72 / D + (2.51 / Re) * ifD )
    sol = root(idarcy, fD0 ** -0.5, args=(ep, D, Re),)
    return  sol.x[0] ** -2

def calculate_fanning(ep, D, Re):
    """Calculate the fanning friction factor"""
    A = (2.457*np.log(((7/Re)**0.9+0.27*ep/D)**-1))**16
    B = (37530/Re)**16
    ff = 2 *((8/Re)**12 + (A+B)**-1.5) ** (1/12)
    return  ff

LMTD Method equations



In [10]:

    
def calculate_heat(water_flowrate_per_tube,
                   water_heat_capacity,
                   water_inlet_temperature,
                   water_outlet_temperature):
    """Calculate the duty"""
    return water_flowrate_per_tube * water_heat_capacity * (water_outlet_temperature - water_inlet_temperature)


def calculate_LMTD(T1a,T1b,T2a,T2b):
    """Calculate LMTD"""
    dT1 = T1a - T1b
    dT2 = T2a - T2b
    return (dT1 - dT2) / np.log(dT1/dT2)

def calculate_heat_via_LMTD( UA, LMTD):
    """Calculate heat exchanger via LMTD"""
    return UA * LMTD
    
def calculate_heat_balance(x0, exhaust_steam_temperature, water_inlet_temperature, water_flowrate_per_tube, water_heat_capacity, UA ):
    """ Heat balance """
    def heat_balance(x, exhaust_steam_temperature, water_inlet_temperature, water_flowrate_per_tube, water_heat_capacity, UA ):
        water_outlet_temperature = x
        LMTD = calculate_LMTD(exhaust_steam_temperature,water_inlet_temperature,exhaust_steam_temperature,water_outlet_temperature)
        Q1 = calculate_heat_via_LMTD( UA, LMTD)
        Q2 = calculate_heat( water_flowrate_per_tube, water_heat_capacity, water_inlet_temperature, water_outlet_temperature)
        return Q1 - Q2
    sol = root(heat_balance, x0, args=(exhaust_steam_temperature, water_inlet_temperature, water_flowrate_per_tube, water_heat_capacity, UA ))
    return  sol.x[0]

Pump parameters



In [11]:

    
def calc_parabola_vertex(x1, y1, x2, y2):
    '''Adapted and modifed to get the unknowns for defining a parabola:
    http://stackoverflow.com/questions/717762/how-to-calculate-the-vertex-of-a-parabola-given-three-points'''
    A = (y1 - y2)/(x1**2 - x2**2)
    B = 0
    C = y1 - A * x1**2 
    return A, B, C

def calc_pump_pressure(x,A,B,C):
        return A*x**2+B*x+C

Pump



In [12]:

    
pump_A, pump_B, pump_C, pump_Kac = [], [], [], []
case_i = 0
for x,y,water_desired_speed in zip(nodeA_w,nodeA_P,v_cases):
    pump_x1 = x - 50
    pump_y1 = y + 10000
    pump_x2 = x
    pump_y2 = y
    pump_Ai,pump_Bi,pump_Ci = calc_parabola_vertex(pump_x1, pump_y1, pump_x2, pump_y2)
    p_calc = calc_pump_pressure(x,pump_Ai,pump_Bi,pump_Ci)
    dp = p_calc - y
    Kac = 2 * dp / (water_density * water_desired_speed**2)
    pump_A.append(pump_Ai)
    pump_B.append(pump_Bi)
    pump_C.append(pump_Ci)
    pump_Kac.append(Kac)
    print("Pump model for case {} is: p = {}*w^2 + ({})*w + ({})".format(case_i,pump_Ai,pump_Bi,pump_Ci))
    case_i += 1









    



Pump model for case 0 is: p = -0.5648760097158674*w^2 + (0)*w + (127117.54750607241)
Pump model for case 1 is: p = -0.5014290728576444*w^2 + (0)*w + (145185.89317053603)
Pump model for case 2 is: p = -0.4289084280506113*w^2 + (0)*w + (135099.81776753164)
Pump model for case 3 is: p = -0.3747844989131249*w^2 + (0)*w + (140178.2403118207)

Preparing cases/scenarios



In [13]:

    
for model in ("simplified","complete"):
    
    # Defining the condenser constructor according to the model type
    if model == "simplified":
        json_condenser = json_condenser_simplified
    else:
        json_condenser = json_condenser_complete

    # Collecting parameters according to the model
    Ntubes_total = tube_bundle[model]["Ntubes_total"]
    Ntubes = tube_bundle[model]["Ntubes"]
    Npipes = tube_bundle[model]["Npipes"]
    print("Ntubes is {number:.{digits}f}".format(number=Ntubes,digits=2))
    print("Npipes is {number:.{digits}f}".format(number=Npipes,digits=2))
    print("Ntubes_total is {number:.{digits}f}".format(number=Ntubes_total,digits=2))
    print("Total area is {number:.{digits}f} m2".format(number=Acold*Ntubes_total, digits=2))

    # This is the file path structure
    filename_structure = "../cases/{case}/scenario_{scenario}/model_{model}/bc_{bc}/{instance}.json"

    case_i = 0
    for water_desired_speed, pump_Ai, pump_Bi, pump_Ci, Tbf, vbf, estimated_steam_mass_flowrate \
        in zip(v_cases,pump_A,pump_B,pump_C,Tbfs,vbfs, kvap_cases):

        print("#####CASE{}######".format(case_i))

        # Selecting required cases
        if water_desired_speed not in required_cases:
            case_i += 1
            print("skiped")
            continue
        
        # Hydraulics
        water_flowrate_per_tube = water_density * Ai * water_desired_speed
        water_flowrate_total = Ntubes_total * water_flowrate_per_tube

        # Calculate head loss of header
        A_header = 0.25 * pi * D_header ** 2
        header_velocity = water_flowrate_total/ (water_density * A_header)
        Re_header = D_header * water_density * header_velocity / water_viscosity
        ff_header = calculate_fanning(ep, D_header, Re_header)
        fD_header = 4 * ff_header
        hL = 0.5 * fD_header * header_velocity ** 2 /  (D_header * g )
        DeltaPOfNozzleInlet = 0.5 * KnozzleInlet * water_density * header_velocity ** 2
        DeltaPOfNozzleOutlet = 0.5 * KnozzleOutlet * water_density * header_velocity ** 2
        DeltaPOfHeader = g * water_density * hL * header_L

        # Calculate head loss of tubes
        Re_tube = Di * water_density * water_desired_speed / water_viscosity
        ff = calculate_fanning(ep, Di, Re_tube)
        g = 9.81
        fD_tube = 4 * ff
        tau = (1/8)*fD_tube*water_density*water_desired_speed**2
        hL = 0.5 * fD_tube * water_desired_speed ** 2 /  (Di * g )
        DeltaP = g * water_density * hL
        pressure_loss_in_tube = DeltaP * L
        DeltaPOfTubeInlet = 0.5 * KtubeInlet * water_density * water_desired_speed ** 2
        DeltaPOfTubeOutlet = 0.5 * KtubeOutlet * water_density * water_desired_speed ** 2

        # Calculating the pressure at the nodes
        pressure_node_D = water_outlet_pressure
        pressure_node_C = pressure_node_D + DeltaPOfHeader + DeltaPOfNozzleOutlet
        pressure_node_B = pressure_node_C + DeltaPOfTubeOutlet + DeltaPOfTubeInlet + pressure_loss_in_tube
        pressure_node_A = pressure_node_B + DeltaPOfNozzleInlet + DeltaPOfHeader

        # Wall conduction
        Reswall = np.log(Do / Di) / (2 * pi * kwall * L)

        # Iterative loop for solving the LMTD with variable properties
        ii = 0
        while ii<10:

            # Internal convection coefficient

            # The properties are calculated at the mean temperature between bulk and internal surface
            int_film_properties_temperature = 0.5*(Tm + Ti)
            int_film_properties_pressure  = Pm

            # Calculating properties
            int_film_properties = _Region1(int_film_properties_temperature,1e-6*int_film_properties_pressure)
            int_film_density = 1/int_film_properties['v']
            int_film_entalphy = 1e3*int_film_properties['h']
            int_film_heat_capacity = 1e3*int_film_properties['cp']
            int_film_conductivity = _ThCond(int_film_density, int_film_properties_temperature)
            int_film_viscosity = _Viscosity(int_film_density, int_film_properties_temperature)
            int_film_prandtl = int_film_heat_capacity * int_film_viscosity / int_film_conductivity

            # Calculating Darcy
            int_film_Re = Di * int_film_density * water_desired_speed / int_film_viscosity
            int_film_fD = calculate_darcy(fD_tube, ep, Di, int_film_Re)

            # Calculates the Nussel dimensionless number using Petukhov correlation modified by Gnielinski. See Incropera 4th Edition [8.63]
            nusselt = (int_film_fD / 8.) * (int_film_Re - 1000.) * int_film_prandtl / ( 1. + 12.7 * (int_film_fD / 8.) ** 0.5 * (int_film_prandtl ** (2 / 3) - 1.))
            #nusselt = 0.023 * Re ** (4/5) * water_prandtl ** 0.4

            hint = nusselt * water_conductivity / Di
            Resint = 1 / (pi * Di * hint * L)


            # External convection coefficient

            # The properties are calculated at the mean temperature between bulk and external surface
            film_properties_temperature = 0.5*(exhaust_steam_temperature + To)
            film_properties_pressure  = exhaust_steam_pressure

            # Calculating properties
            film_properties = _Region1(film_properties_temperature,1e-6*film_properties_pressure)
            film_density = 1/film_properties['v']
            film_entalphy = 1e3*film_properties['h']
            film_heat_capacity = 1e3*film_properties['cp']
            film_conductivity = _ThCond(film_density, film_properties_temperature)
            film_viscosity = _Viscosity(film_density, film_properties_temperature)
            film_prandtl = film_heat_capacity * film_viscosity / film_conductivity

            num = (g * film_density * (film_density - exhaust_steam_density) * film_conductivity ** 3. * exhaust_vaporization_heat)
            den = film_viscosity * (exhaust_steam_temperature - To) * Do
            hext = 0.729 * (num / den) ** 0.25
            Resext = 1 / (pi * Do * hext * L)

            # Effect of vertical position
            tube_index = np.linspace(1,Ntubes,num=Ntubes)
            if model == "simplified":
                fNtub = 0*tube_index + tubes_per_width**(-1/6)
                fNtub_soft = fNtub
                hext_list = hext*fNtub
                hext_list_soft = np.mean(hext*fNtub_soft)
            else:
                fNtub = (tube_index ** (5/6) - (tube_index - 1) ** (5/6))
                fNtub_soft = fNtub
                #fNtub_soft = (tube_index ** 0.95 - (tube_index - 1) ** 0.95)
                hext_list = hext*fNtub
                hext_list_soft = hext*fNtub_soft
            hext_list = np.tile(hext_list,(Nelements,1))
            hext_list_soft = np.tile(hext_list_soft,(Nelements,1))

            
            # Simulate heat exchanger

            # Calculate Overall he
            Restotal = Resext + Resint + Reswall
            UA = 1/Restotal
            Ucold = UA / Acold
            Uhot = UA / Ahot

            # Initial Estimate of water_outlet_temperature
            x0 = water_inlet_temperature + 10

            # Calculate water_outlet_temperature
            water_outlet_temperature = calculate_heat_balance(x0, exhaust_steam_temperature, water_inlet_temperature, water_flowrate_per_tube, water_heat_capacity, UA )

            # Calculate LMTD
            LMTD = calculate_LMTD(exhaust_steam_temperature,water_inlet_temperature,exhaust_steam_temperature,water_outlet_temperature)

            Q = calculate_heat_via_LMTD( UA, LMTD)

            Q_at_inlet = (exhaust_steam_temperature - water_inlet_temperature) / Restotal
            temperature_internal_surface_at_inlet = water_inlet_temperature + Q_at_inlet*Resint
            temperature_external_surface_at_intlet = exhaust_steam_temperature - Q_at_inlet*Resext

            Q_at_outlet = (exhaust_steam_temperature - water_outlet_temperature) / Restotal
            temperature_internal_surface_at_outlet = water_outlet_temperature + Q_at_outlet*Resint
            temperature_external_surface_at_outlet = exhaust_steam_temperature - Q_at_outlet*Resext

            # Confirming the surface temperature estimates
            TmNew = 0.5 * (water_inlet_temperature + water_outlet_temperature)
            TiNew = 0.5 * (temperature_internal_surface_at_inlet + temperature_internal_surface_at_outlet)
            ToNew = 0.5 * (temperature_external_surface_at_intlet + temperature_external_surface_at_outlet)

            TmOld, TiOld, ToOld = Tm, Ti, To
            Tm, Ti, To = TmNew, TiNew, ToNew

            if (ToNew - ToOld)**2 < 1e-6:
                print("Conderged!")
                break
                
            ii+=1
            
        # Calculating condensation
        A_total = Acold*n_tubes
        Q_total = Q*n_tubes
        steam_mass_flowrate =  Q_total/exhaust_vaporization_heat

        print("The water velocity is {number:.{digits}f} m/s".format(number=water_desired_speed, digits=3))
        print("The water flow rate is {number:.{digits}f} kg/s per tube".format(number=water_flowrate_per_tube, digits=5))
        print("The total water flow rate is {number:.{digits}f} kg/s".format(number=water_flowrate_total, digits=5))
        print("Estimated internal surface temperature at inlet is {number:.{digits}f} K".format(number=Ti0,digits=2))
        print("Estimated external surface temperature at inlet is {number:.{digits}f} K".format(number=To0,digits=2))
        print("Estimated internal surface temperature at outlet is {number:.{digits}f} K".format(number=Ti1,digits=2))
        print("Estimated external surface temperature at outlet is {number:.{digits}f} K".format(number=To1,digits=2))
        print("The total water flowrate is {number:.{digits}f} ton/h".format(number=water_flowrate_total*3.6, digits=2))
        print("The Header Reynolds number is {number:.{digits}f}.".format(number=Re_header, digits=0))
        print("The Header velocity is {number:.{digits}f} m/s.".format(number=header_velocity, digits=0))
        print("The header diameter is {number:.{digits}f} in".format(number=D_header/0.0254, digits=2))
        print("The Pressure loss of the inlet nozzle is {number:.{digits}f} Pa.".format(number=DeltaPOfNozzleInlet, digits=0))
        print("The Pressure loss of the outlet nozzle is {number:.{digits}f} Pa.".format(number=DeltaPOfNozzleOutlet, digits=0))
        print("The Pressure loss of the header is {number:.{digits}f} Pa.".format(number=DeltaPOfHeader, digits=0))
        print("The Reynolds number is {number:.{digits}f}.".format(number=Re_tube, digits=0))
        print("The calculated friction factor is: {number:.{digits}f}".format(number=fD_tube, digits=4))
        print("The wall shear stress is: {number:.{digits}f} Pa".format(number=tau, digits=4))
        print("The head loss is equal to {number:.{digits}f} m/m".format(number=hL, digits=4))
        print("The pressure loss is equal to {number:.{digits}f} Pa/m".format(number=DeltaP, digits=1))
        print("The pressure loss in tube is equal to {number:.{digits}f} Pa".format(number=pressure_loss_in_tube, digits=1))
        print("The Pressure loss of the inlet tube is {number:.{digits}f} Pa.".format(number=DeltaPOfTubeInlet, digits=0))
        print("The Pressure loss of the outlet tube is {number:.{digits}f} Pa.".format(number=DeltaPOfTubeOutlet, digits=0))
        print("The Pressure at node A is {number:.{digits}f} Pa.".format(number=pressure_node_A, digits=0))
        print("The Pressure at node B is {number:.{digits}f} Pa.".format(number=pressure_node_B, digits=0))
        print("The Pressure at node C is {number:.{digits}f} Pa.".format(number=pressure_node_C, digits=0))
        print("The Pressure at node D is {number:.{digits}f} Pa.".format(number=pressure_node_D, digits=0))
        print("The calculated internal convection coefficient is {number:.{digits}f} W/(K*m2)".format(number=hint,digits=2))
        print("The calculated internal resistance x length is {number:.{digits}f} K/W".format(number=Resint,digits=6))
        print("The calculated external convection coefficient is {number:.{digits}f} W/(K*m2)".format(number=hext,digits=2))
        print("The calculated external resistance x length is {number:.{digits}f} K/W".format(number=Resext,digits=6))
        print("The overall resistance is {number:.{digits}f} K/W".format(number=Restotal,digits=6))
        print("U*A per tube is {number:.{digits}f} W/K".format(number=UA,digits=2))
        print("The tube for cold side is {number:.{digits}f} m2".format(number=Acold,digits=2))
        print("The tube for hot side is {number:.{digits}f} m2".format(number=Ahot,digits=2))
        print("The overall heat exchange coefficient U for cold side is {number:.{digits}f} W/(K*m2)".format(number=Ucold,digits=2))
        print("The overall heat exchange coefficient U for hot side is {number:.{digits}f} W/(K*m2)".format(number=Uhot,digits=2))
        print("LMTD is {number:.{digits}f} K".format(number=LMTD,digits=2))
        print("The cooling water outlet temperature is {number:.{digits}f} K".format(number=water_outlet_temperature,digits=2))
        print("The cooling water inlet temperature is {number:.{digits}f} K".format(number=water_inlet_temperature,digits=2))
        print("The heat rate per tube is {number:.{digits}f} W".format(number=Q,digits=2))
        print("The heat rate per tube length is {number:.{digits}f} W".format(number=Q/L,digits=2))
        print("Heat rate per tube length at inlet is {0} W/m".format(Q_at_inlet/L))
        print("Heat rate per tube length at outlet is {0} W/m".format(Q_at_outlet/L))
        print("Pipe temperature at internal surface for the inlet is {number:.{digits}f} K".format(number=temperature_internal_surface_at_inlet,digits=2))
        print("Pipe temperature at external surface for the inlet is {number:.{digits}f} K".format(number=temperature_external_surface_at_intlet,digits=2))
        print("Pipe temperature at internal surface for the outlet is {number:.{digits}f} K".format(number=temperature_internal_surface_at_outlet,digits=2))
        print("Pipe temperature at external surface for the outlet is {number:.{digits}f} K".format(number=temperature_external_surface_at_outlet,digits=2))
        print("The total transfer area (cold side) is {number:.{digits}f} m2".format(number=A_total, digits=2))
        print("The total steam flowrate is {number:.{digits}f} ton/h".format(number=steam_mass_flowrate*3.6, digits=2))
        print("The total heat load of condensation is {number:.{digits}f} MW".format(number=1e-6*Q_total, digits=2))
        print("The calculated wall resistance x length is {number:.{digits}f} K/W".format(number=Reswall,digits=6))

        # Preparing JSON fracments
        
        # nodes
        json_nodeA = json_node_inlet(pressure_node_A,water_inlet_temperature, water_flowrate_total)
        json_nodeAw = json_node_wfix(pressure_node_A,water_inlet_temperature, water_flowrate_total)
        json_nodeA_pumped = json_node_pumped(pressure_node_A,water_inlet_temperature, water_flowrate_total, pump_Ai, pump_Bi, pump_Ci)
        json_nodeA_river1 = json_node_river(pressure_node_A,water_inlet_temperature, water_flowrate_total, pump_Ai, pump_Bi, pump_Ci, 13.8+273.15, 6.2, -121.0, "Fixed")
        json_nodeA_river2 = json_node_river(pressure_node_A,water_inlet_temperature, water_flowrate_total, pump_Ai, pump_Bi, pump_Ci, 13.8+273.15, 6.2, -121.0, "PreVariable")
        json_nodeB = json_node_wfix(pressure_node_B,water_inlet_temperature, 0.0)
        json_nodeC = json_node_wfix(pressure_node_C,water_inlet_temperature, 0.0)
        json_nodeD = json_node_outlet(pressure_node_D, water_inlet_temperature, water_flowrate_total)

        # pipes
        json_pipe1 = json_pipe("node_A", "node_B", Nelements, D_header, header_L, ep, 0.0, KnozzleInlet, Re_header, header_velocity, water_flowrate_total, water_inlet_temperature, pressure_node_A, pressure_node_B, fD_header,)
        json_pipe2 = json_pipe("node_C", "node_D", Nelements, D_header, header_L, ep, KnozzleOutlet, 0.0, Re_header, header_velocity, water_flowrate_total, water_outlet_temperature, pressure_node_C, pressure_node_D, fD_header,)

        # condenser
        cond_json_soft = json_condenser("node_B", "node_C", Nelements, Ntubes, Npipes,\
                       Di, Do, kwall, L, ep,\
                       water_flowrate_per_tube, water_inlet_temperature, water_outlet_temperature,\
                       pressure_node_B, pressure_node_C,\
                       temperature_external_surface_at_intlet, temperature_external_surface_at_outlet,\
                       temperature_internal_surface_at_inlet, temperature_internal_surface_at_outlet,\
                       Re_tube, water_desired_speed,\
                       fD_tube, hint, hext_list_soft,\
                       external_volume,\
                       exhaust_steam_pressure, exhaust_steam_pressure + 20000, exhaust_steam_pressure + 5000,\
                       estimated_steam_mass_flowrate,\
                       fNtub_soft,\
                       exhaust_steam_pressure, exhaust_steam_temperature,\
                       rhomf, Ccell, Csubstrate, lagt_infinite, 'Fixed', vbf, Tbf, Q_at_inlet, Q_at_outlet )

        cond_json = json_condenser("node_B", "node_C", Nelements, Ntubes, Npipes,\
                       Di, Do, kwall, L, ep,\
                       water_flowrate_per_tube, water_inlet_temperature, water_outlet_temperature,\
                       pressure_node_B, pressure_node_C,\
                       temperature_external_surface_at_intlet, temperature_external_surface_at_outlet,\
                       temperature_internal_surface_at_inlet, temperature_internal_surface_at_outlet,\
                       Re_tube, water_desired_speed,\
                       fD_tube, hint, hext_list,\
                       external_volume,\
                       exhaust_steam_pressure, exhaust_steam_pressure + 20000, exhaust_steam_pressure + 5000,\
                       estimated_steam_mass_flowrate,\
                       fNtub,\
                       exhaust_steam_pressure,\
                       exhaust_steam_temperature,\
                       rhomf, Ccell, Csubstrate, lagt, 'PreVariable', vbf, Tbf, Q_at_inlet, Q_at_outlet )

        cond_json_clean = json_condenser("node_B", "node_C", Nelements, Ntubes, Npipes,\
                       Di, Do, kwall, L, ep,\
                       water_flowrate_per_tube, water_inlet_temperature, water_outlet_temperature,\
                       pressure_node_B, pressure_node_C,\
                       temperature_external_surface_at_intlet, temperature_external_surface_at_outlet,\
                       temperature_internal_surface_at_inlet, temperature_internal_surface_at_outlet,\
                       Re_tube, water_desired_speed,\
                       fD_tube, hint, hext_list,\
                       external_volume,\
                       exhaust_steam_pressure, exhaust_steam_pressure + 20000, exhaust_steam_pressure + 5000,\
                       estimated_steam_mass_flowrate,\
                       fNtub,\
                       exhaust_steam_pressure,\
                       exhaust_steam_temperature,\
                       rhomf, Ccell, Csubstrate, lagt_infinite, 'PreVariable', vbf, Tbf, Q_at_inlet, Q_at_outlet )

        filenames={}
        filenames["A1"] = filename_structure.format(case=case_i,scenario="shortterm", model=model, bc="pfixed", instance="preparation")
        filenames["A2"] = filename_structure.format(case=case_i,scenario="shortterm", model=model, bc="wfixed", instance="preparation")
        filenames["A3"] = filename_structure.format(case=case_i,scenario="shortterm", model=model, bc="pump", instance="preparation")
        
        filenames["B1"] = filename_structure.format(case=case_i,scenario="shortterm", model=model, bc="pfixed", instance="simulation")
        filenames["B2"] = filename_structure.format(case=case_i,scenario="shortterm", model=model, bc="wfixed", instance="simulation")
        filenames["B3"] = filename_structure.format(case=case_i,scenario="shortterm", model=model, bc="pump", instance="simulation")
        
        filenames["C3"] = filename_structure.format(case=case_i,scenario="shortterm", model=model, bc="pump", instance="simulationwithoutbiofilm")
        
        filenames["A4"] = filename_structure.format(case=case_i,scenario="longterm", model=model, bc="pump", instance="preparation")
        filenames["B4"] = filename_structure.format(case=case_i,scenario="longterm", model=model, bc="pump", instance="simulation")
        filenames["C4"] = filename_structure.format(case=case_i,scenario="longterm", model=model, bc="pump", instance="simulationwithoutbiofilm")
        
        simulation_parameters={}
        simulation_parameters["A1"] = simulation_parameters_default["shortterm_preparation"]
        simulation_parameters["A2"] = simulation_parameters_default["shortterm_preparation"]
        simulation_parameters["A3"] = simulation_parameters_default["shortterm_preparation"]
        simulation_parameters["B1"] = simulation_parameters_default["shortterm_simulation"]
        simulation_parameters["B2"] = simulation_parameters_default["shortterm_simulation"]
        simulation_parameters["B3"] = simulation_parameters_default["shortterm_simulation"]
        simulation_parameters["C3"] = simulation_parameters_default["shortterm_simulation"]

        simulation_parameters["A4"] = simulation_parameters_default["longterm_preparation"]
        simulation_parameters["B4"] = simulation_parameters_default["longterm_simulation"]
        simulation_parameters["C4"] = simulation_parameters_default["longterm_simulation"]

        json_data_cases={}
        json_data_cases["A1"] = { "class": "Network",
                                 "submodels": { "node_A": json_nodeA, "condenser": cond_json_soft,} }

        json_data_cases["A2"] = { "class": "Network",
                                 "submodels": { "node_A": json_nodeAw, "condenser": cond_json_soft, } }

        json_data_cases["A3"] = { "class": "Network",
                                 "submodels": { "node_A": json_nodeA_pumped, "condenser": cond_json_soft, } }


        json_data_cases["B1"] = { "class": "Network",
                                 "submodels": { "node_A": json_nodeA, "condenser": cond_json,  } }

        json_data_cases["B2"] = { "class": "Network",
                                 "submodels": { "node_A": json_nodeAw, "condenser": cond_json, } }

        json_data_cases["B3"] = { "class": "Network",
                                 "submodels": { "node_A": json_nodeA_pumped, "condenser": cond_json, } }


        json_data_cases["C3"] = { "class": "Network",
                                 "submodels": { "node_A": json_nodeA_pumped, "condenser": cond_json_clean, } }

        
        json_data_cases["A4"] = { "class": "Network",
                                 "submodels": { "node_A": json_nodeA_river1, "condenser": cond_json_soft,} }

        json_data_cases["B4"] = { "class": "Network",
                                 "submodels": { "node_A": json_nodeA_river2, "condenser": cond_json, } }

        json_data_cases["C4"] = { "class": "Network",
                                 "submodels": {"node_A": json_nodeA_river2, "condenser": cond_json_clean,} }

        for key, value in json_data_cases.items():
            filename_i = filenames[key]
            name_i="case_{}_{}".format(case_i, key)
            print(name_i)
            value["name"] = name_i
            value["simulation_parameters"] = simulation_parameters[key]
            value["submodels"]["node_B"] = json_nodeB
            value["submodels"]["node_C"] = json_nodeC
            value["submodels"]["node_D"] = json_nodeD
            value["submodels"]["pipe_1"] = json_pipe1
            value["submodels"]["pipe_2"] = json_pipe2
            # Creating directory
            directory = os.path.dirname(filename_i)
            if not os.path.exists(directory):
                os.makedirs(directory)
            with open(filename_i, 'w') as outfile:
                json.dump(value, outfile, indent=4)

        case_i += 1









    



Ntubes is 1.00
Npipes is 1156.00
Ntubes_total is 1156.00
Total area is 261.51 m2
#####CASE0######
skiped
#####CASE1######
Conderged!
The water velocity is 1.000 m/s
The water flow rate is 0.19416 kg/s per tube
The total water flow rate is 224.44565 kg/s
Estimated internal surface temperature at inlet is 312.00 K
Estimated external surface temperature at inlet is 314.00 K
Estimated internal surface temperature at outlet is 312.00 K
Estimated external surface temperature at outlet is 314.00 K
The total water flowrate is 808.00 ton/h
The Header Reynolds number is 941310.
The Header velocity is 2 m/s.
The header diameter is 14.00 in
The Pressure loss of the inlet nozzle is 2819 Pa.
The Pressure loss of the outlet nozzle is 1794 Pa.
The Pressure loss of the header is 5024 Pa.
The Reynolds number is 18385.
The calculated friction factor is: 0.0320
The wall shear stress is: 3.9894 Pa
The head loss is equal to 0.1036 m/m
The pressure loss is equal to 1013.2 Pa/m
The pressure loss in tube is equal to 4632.3 Pa
The Pressure loss of the inlet tube is 224 Pa.
The Pressure loss of the outlet tube is 224 Pa.
The Pressure at node A is 119742 Pa.
The Pressure at node B is 111899 Pa.
The Pressure at node C is 106818 Pa.
The Pressure at node D is 100000 Pa.
The calculated internal convection coefficient is 5838.90 W/(K*m2)
The calculated internal resistance x length is 0.000757 K/W
The calculated external convection coefficient is 12784.87 W/(K*m2)
The calculated external resistance x length is 0.000286 K/W
The overall resistance is 0.001175 K/W
U*A per tube is 850.80 W/K
The tube for cold side is 0.23 m2
The tube for hot side is 0.27 m2
The overall heat exchange coefficient U for cold side is 3760.90 W/(K*m2)
The overall heat exchange coefficient U for hot side is 3109.40 W/(K*m2)
LMTD is 19.83 K
The cooling water outlet temperature is 307.74 K
The cooling water inlet temperature is 286.95 K
The heat rate per tube is 16873.15 W
The heat rate per tube length is 3690.54 W
Heat rate per tube length at inlet is 5956.271696958621 W/m
Heat rate per tube length at outlet is 2088.381146624725 W/m
Pipe temperature at internal surface for the inlet is 307.57 K
Pipe temperature at external surface for the inlet is 311.17 K
Pipe temperature at internal surface for the outlet is 314.96 K
Pipe temperature at external surface for the outlet is 316.23 K
The total transfer area (cold side) is 261.51 m2
The total steam flowrate is 29.35 ton/h
The total heat load of condensation is 19.51 MW
The calculated wall resistance x length is 0.000132 K/W
case_1_A1
case_1_A2
case_1_A3
case_1_B1
case_1_B2
case_1_B3
case_1_C3
case_1_A4
case_1_B4
case_1_C4
#####CASE2######
skiped
#####CASE3######
skiped
Ntubes is 34.00
Npipes is 34.00
Ntubes_total is 1156.00
Total area is 261.51 m2
#####CASE0######
skiped
#####CASE1######
Conderged!
The water velocity is 1.000 m/s
The water flow rate is 0.19416 kg/s per tube
The total water flow rate is 224.44565 kg/s
Estimated internal surface temperature at inlet is 312.00 K
Estimated external surface temperature at inlet is 314.00 K
Estimated internal surface temperature at outlet is 312.00 K
Estimated external surface temperature at outlet is 314.00 K
The total water flowrate is 808.00 ton/h
The Header Reynolds number is 941310.
The Header velocity is 2 m/s.
The header diameter is 14.00 in
The Pressure loss of the inlet nozzle is 2819 Pa.
The Pressure loss of the outlet nozzle is 1794 Pa.
The Pressure loss of the header is 5024 Pa.
The Reynolds number is 18385.
The calculated friction factor is: 0.0320
The wall shear stress is: 3.9894 Pa
The head loss is equal to 0.1036 m/m
The pressure loss is equal to 1013.2 Pa/m
The pressure loss in tube is equal to 4632.3 Pa
The Pressure loss of the inlet tube is 224 Pa.
The Pressure loss of the outlet tube is 224 Pa.
The Pressure at node A is 119742 Pa.
The Pressure at node B is 111899 Pa.
The Pressure at node C is 106818 Pa.
The Pressure at node D is 100000 Pa.
The calculated internal convection coefficient is 5838.87 W/(K*m2)
The calculated internal resistance x length is 0.000757 K/W
The calculated external convection coefficient is 12784.47 W/(K*m2)
The calculated external resistance x length is 0.000286 K/W
The overall resistance is 0.001175 K/W
U*A per tube is 850.79 W/K
The tube for cold side is 0.23 m2
The tube for hot side is 0.27 m2
The overall heat exchange coefficient U for cold side is 3760.86 W/(K*m2)
The overall heat exchange coefficient U for hot side is 3109.37 W/(K*m2)
LMTD is 19.83 K
The cooling water outlet temperature is 307.73 K
The cooling water inlet temperature is 286.95 K
The heat rate per tube is 16873.05 W
The heat rate per tube length is 3690.52 W
Heat rate per tube length at inlet is 5956.209215838814 W/m
Heat rate per tube length at outlet is 2088.3821992655403 W/m
Pipe temperature at internal surface for the inlet is 307.57 K
Pipe temperature at external surface for the inlet is 311.17 K
Pipe temperature at internal surface for the outlet is 314.96 K
Pipe temperature at external surface for the outlet is 316.23 K
The total transfer area (cold side) is 261.51 m2
The total steam flowrate is 29.35 ton/h
The total heat load of condensation is 19.51 MW
The calculated wall resistance x length is 0.000132 K/W
case_1_A1
case_1_A2
case_1_A3
case_1_B1
case_1_B2
case_1_B3
case_1_C3
case_1_A4
case_1_B4
case_1_C4
#####CASE2######
skiped
#####CASE3######
skiped