search

JungHulk / Hulk-Engineering

1 stars 0 forks source link

Vent line heat transfer #2

Open JungHulk opened 2 years ago

JungHulk commented 2 years ago 
#!/usr/bin/env python
# coding: utf-8

# # Calculation of Pipe Heat transfer
# ## Consider pipe cooling and LNG evaporation
# 

# In[1]:

import os; os.environ['RPPREFIX'] = 'C:\\Program Files (x86)\\REFPROP'
from ctREFPROP.ctREFPROP import REFPROPFunctionLibrary
RP = REFPROPFunctionLibrary(os.environ['RPPREFIX'])
RP.SETPATHdll(os.environ['RPPREFIX'])

# In[2]:

from win32com.client import Dispatch
import os
import time; xxstart = time.time()

# In[3]:

import math
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image

openfiledir = os.getcwd()  # 현재 폴더 위치 

# In[4]:

MASS_BASE_SI = RP.GETENUMdll(0, "MASS BASE SI").iEnum   # output: 21

# In[8]:

lng_density = 430                     # [kg/m3]
lng_latent_heat = 510000              # [J/kg]

time_release = 1               # [s]
time_total = 100                # [s]
delta_time = 0.1                # [s]

time_number = int(time_total/delta_time)         # [#].
times = np.linspace(0,time_total, time_number)

inner_diameter = 0.1143                 # [m]  50A
pipe_thickness = 0.003048                 # [m]  Sch. 10
# inner_diameter = 0.24448                 # [m]  50A
# pipe_thickness = 0.008687                 # [m]  Sch. 10
outer_diameter = inner_diameter + 2*pipe_thickness                                          # [m]
cross_area = 0.25*math.pi*inner_diameter**2                                                 # [m2]

# 2.53 kg, 36 초간 release
flow_rate = np.zeros(time_number) 
masstot = np.zeros(time_number) 
enthaltot = np.zeros(time_number) 
# flow_rate[:int(time_release/delta_time)] = 4535/lng_density/time_release/3600               # [m3/s]
flow_rate[:int(time_release/delta_time)] = 34790/lng_density/time_release/3600               # [m3/s]
# flow_rate[:int(time_total/delta_time)] = 2.53/lng_density/time_release                    # [m3/s]
flow_velocity = flow_rate[0]/cross_area
distance_progress = flow_velocity * delta_time

# In[9]:

delta_length = distance_progress                   # [m]
node_number = 2000    #2000

line_length =  delta_length * node_number          # [m]
distance_node = np.linspace(0, line_length, node_number)

index_end_variable = 0

# In[10]:

ambient_temperature = 30 + 273.15                  # [K]
ambient_temp_node = ambient_temperature*np.ones(node_number)

# In[11]:

# In[12]:

pipe_volume_unit = 0.25*math.pi*(outer_diameter**2 - inner_diameter**2)*delta_length     # [m3]

inner_surface_tot = math.pi*inner_diameter*line_length                                   # [m2]
area_unit_in = math.pi*inner_diameter*delta_length                                       # [m2]

fluid_vol_unit = cross_area * delta_length
inner_surface_tot, area_unit_in, cross_area

# In[13]:

# In[14]:

# pipe material: SUS316L
pipe_thermal_conductivity = 16                     # [W/mk]
pipe_density = 8030                                # [kg/m3]
pipe_specific_heat = 3.1121*ambient_temperature - 0.0051*ambient_temperature**2                          # [J/kg K]

# In[15]:

pipe_heat_capa_unit = pipe_specific_heat*pipe_density*pipe_volume_unit                   # [J/K]
pipe_heat_capa_unit

# In[16]:

insulation_thickness = 0.03                                                                               # [m]
insulation_outer_diameter = outer_diameter + 2 * pipe_thickness                                           # [m]

insulation_volume_unit = 0.25*math.pi*(insulation_outer_diameter**2 - outer_diameter**2)*delta_length     # [m3]
insulation_volume_unit, pipe_thickness

area_unit_ex = math.pi*insulation_outer_diameter*delta_length 

# In[17]:

# Insulation material: Rock Wool
insulation_density = 55                                                                                   # [kg/m3]
insulation_specific_heat = 840                                                                            # [J/kg K]
insulation_thermal_conductivity = 0.03                                                                    # [W/mK]

insulation_heat_capa_unit = insulation_specific_heat*insulation_density*insulation_volume_unit            # [J/K]
insulation_heat_capa_unit

# In[18]:

hc_in= 800                                                                      # [W/m2K]
hc_ex = 8                                                                      # [W/m2K]
pipe_heat_coeff = pipe_thermal_conductivity/pipe_thickness                                # [W/m2K]
insulation_heat_coeff = insulation_thermal_conductivity/insulation_thickness              # [W/m2K]

# overall_heat_coeff = 1/(1/inner_heat_coeff + 1/pipe_heat_coeff + 1/insulation_heat_coeff +1/outer_heat_coeff)
overall_heat_coeff = 1/(1/hc_in + 1/pipe_heat_coeff +1/hc_ex)
overall_heat_coeff

# In[19]:

################### fluid properties #########################
fluid_temp_ini = -120 + 273.15                                                      # [K]
pressure = (16.2 + 1)* 101300                                                               # [Pa]

# In[22]:

fluid = "methane"
znit = [1.]
flag_node = np.zeros(node_number)

lng_density = RP.REFPROPdll(fluid, "PT", "D", MASS_BASE_SI, 1, 0, pressure, fluid_temp_ini, znit).Output[0]         # [kg/m3]
lng_specific_heat = RP.REFPROPdll(fluid, "PT", "Cp", MASS_BASE_SI, 1, 0, pressure, fluid_temp_ini, znit).Output[0]  # [J/kg K]

ng_density = RP.REFPROPdll(fluid, "PQ", "D", MASS_BASE_SI, 1, 0, pressure, 1, znit).Output[0]         # [kg/m3]
ng_specific_heat = RP.REFPROPdll(fluid, "PQ", "Cp", MASS_BASE_SI, 1, 0, pressure, 1, znit).Output[0]  # [J/kg K]

specific_heat_ini = lng_specific_heat
density_ini = lng_density 

ng_density, ng_specific_heat, lng_density, lng_specific_heat
index_node = np.linspace(0,node_number-1,node_number)
index_node_incre= np.zeros((node_number, time_number))

# In[23]:

index_history = np.zeros(time_number)
heatflow_ex = np.zeros(time_number)
heatflow_in = np.zeros(time_number)

passing_vapor = np.zeros(time_number)
passing_liquid = np.zeros(time_number)

pipe_temp_node = np.zeros((node_number, time_number))
pipe_specific_heat_node = np.zeros((node_number, time_number))
pipe_heat_capa_node = np.zeros((node_number, time_number))

insulation_temp_node = np.zeros((node_number, time_number))
fluid_temp_node = np.zeros((node_number, time_number))
fluid_quality_node = np.zeros((node_number, time_number))
fluid_density_node = np.zeros((node_number, time_number))
fluid_mass_node = np.zeros((node_number, time_number))
fluid_mass_node2 = np.zeros((node_number, time_number))
fluid_specific_heat_node = np.zeros((node_number, time_number))
fluid_evaporated_vol = np.zeros((node_number, time_number))

heat_flow_ex = np.zeros((node_number, time_number))
heat_flow_in = np.zeros((node_number, time_number))

heat_flux = np.zeros((node_number, time_number))
pipe_node_new = np.zeros((node_number, time_number))

## Initialization
pipe_temp_node[:,0] = ambient_temperature                                                    # [K]
fluid_temp_node[:,:] = ambient_temperature  
fluid_quality_node[:,0] = 0  
fluid_density_node[:,0] = 0
fluid_specific_heat_node[:,0] = lng_specific_heat

# insulation_temp_node[:] = ambient_temperature                                              # [K]

# In[24]:

def Heat_trasnfer_ex(pipe_temp_node, ambient_temp_node, hc_ex, area_unit_ex, delta_time):
    heat_flow_ex = hc_ex*area_unit_ex*(ambient_temp_node - pipe_temp_node) * delta_time

    heat_flow_ex_sum = np.sum(heat_flow_ex)

    results = [heat_flow_ex, heat_flow_ex_sum]
    return results

# In[25]:

def Heat_trasnfer_in(fluid_temp_node, pipe_temp_node, fluid_quality, hc_in, area_unit_in, delta_time, node_number):

    hc_in_node = np.zeros(node_number)
    hc_in_node[:] = hc_in

    for i in range(0, node_number):

        if pipe_temp_node[i] - fluid_temp_node[i] <5:
            hc_in_node[i] = 1500

        if fluid_quality[i] == 1:
            hc_in_node[i] = 20

    heat_flow_in = hc_in_node*area_unit_in*(pipe_temp_node - fluid_temp_node) * delta_time

    heat_flow_in_sum = np.sum(heat_flow_in)

    results = [heat_flow_in, heat_flow_in_sum]
    return results

def Update_pipe_pro(pipe_density, pipe_temp_node, pipe_volume_unit, node_number):

    pipe_specific_heat_new = np.zeros(node_number)
    pipe_heat_capa_new = np.zeros(node_number)

    pipe_specific_heat_new = 3.1121*pipe_temp_node - 0.0051*pipe_temp_node**2
    pipe_heat_capa_new = pipe_specific_heat_new*pipe_density*pipe_volume_unit

    results = [pipe_specific_heat_new, pipe_heat_capa_new]

    return results

# In[31]:

def Upgrade_pipe_temp(heat_flow_ex, heat_flow_in, pipe_specific_heat_node, pipe_temp_node, node_number):
    temp_delt = np.zeros(node_number)

    temp_delt = (heat_flow_ex - heat_flow_in)/pipe_specific_heat_node
    pipe_temp_node = pipe_temp_node + temp_delt

    return pipe_temp_node

# In[36]:

def Upgrade_fluid_properties(heat_flow_in, fluid_mass_node, fluid_temp_node, fluid_quality_node, fluid_density_node, fluid_specific_heat_node, fluid_evaporated_vol, pressure, fluid_vol_unit, flag_node, lng_density, ng_density, lng_specific_heat, ng_specific_heat):
# calculating the mass based fluid properties
    fluid_heat_capa = np.zeros(node_number)
    fluid_temp_new = fluid_temp_node
    fluid_specific_heat_new = np.zeros(node_number)
    fluid_density_new = np.zeros(node_number)
    fluid_mass_new  = np.zeros(node_number)
    fluid_quality_new = np.zeros(node_number)
    fluid_evaporated_vol_new = np.zeros(node_number)

    for i in range(0, node_number):
        mass_evap = 0
        # if i == 175:
        #     aa = 0
        if flag_node[i] >0:

            if fluid_quality_node[i] == 0:      # Liquid case

                mass_evap = heat_flow_in[i]/lng_latent_heat
                mass_tot = fluid_mass_node[i]

                fluid_evaporated_vol_new[i] = fluid_evaporated_vol[i] + mass_evap/ng_density        # [m3]

                if mass_tot ==0:
                    fluid_quality_new[i] = 0
                else:
                    fluid_quality_new[i] = mass_evap/mass_tot
                fluid_temp_new[i] = fluid_temp_node[i]

                fluid_specific_heat_new[i] = ng_specific_heat*fluid_quality_new[i] + lng_specific_heat*(1-fluid_quality_new[i])
                fluid_density_new[i] = fluid_mass_node[i]/fluid_vol_unit
                fluid_mass_new[i] = fluid_mass_node[i]

            elif fluid_quality_node[i] >= 1:      # Vapor case

                fluid_quality_node[i] = 1         
                fluid_heat_capa[i] = fluid_density_node[i]*fluid_specific_heat_node[i]*fluid_vol_unit

                delta_temp = heat_flow_in[i]/fluid_heat_capa[i]
                fluid_evaporated_vol_new[i] = fluid_evaporated_vol[i]

                fluid_temp_new[i] = fluid_temp_node[i] + delta_temp
                fluid_specific_heat_new[i] = ng_specific_heat
                fluid_density_new[i] = fluid_mass_node[i]/fluid_vol_unit
                fluid_quality_new[i] = fluid_quality_node[i]
                fluid_mass_new[i] = fluid_mass_node[i]

            elif fluid_quality_node[i] < 1 and  fluid_quality_node[i] > 0 :                               # Two phase case

                mass_evap = heat_flow_in[i]/lng_latent_heat
                mass_tot = fluid_mass_node[i]   

                fluid_evaporated_vol_new[i] = fluid_evaporated_vol[i] + mass_evap/ng_density
                if mass_tot ==0:
                    fluid_quality_delta = 0
                else:
                    fluid_quality_delta = mass_evap/mass_tot
                fluid_quality_new[i] = fluid_quality_node[i] + fluid_quality_delta

                if fluid_quality_new[i] > 1:
                    fluid_quality_new[i] = 1

                fluid_temp_new[i] = fluid_temp_node[i]            
                fluid_specific_heat_new[i] = ng_specific_heat*fluid_quality_new[i] + lng_specific_heat*(1-fluid_quality_new[i])
                fluid_density_new[i] = fluid_mass_node[i]/fluid_vol_unit          
                fluid_mass_new[i] = fluid_mass_node[i]
        else:
            fluid_temp_new[i] = fluid_temp_node[i]
            fluid_specific_heat_new[i] = fluid_specific_heat_node[i]
            fluid_density_new[i] = fluid_density_node[i]
            fluid_quality_new[i] = fluid_quality_node[i]
            fluid_evaporated_vol_new[i] = fluid_evaporated_vol[i]
            fluid_mass_new[i] = fluid_mass_node[i]         

    results = [fluid_temp_new, fluid_mass_new, fluid_specific_heat_new, fluid_density_new, fluid_quality_new, fluid_evaporated_vol_new] 

    return results

def Vap_expansion(fluid_quality_node, fluid_mass_node, fluid_evaporated_vol, node_number, ng_density, lng_density, fluid_vol_unit,  index_node):

    index_node_new = np.zeros(node_number)
    fluid_evaporated_vol_new = fluid_evaporated_vol
    index_node_incre = np.zeros(node_number)  
    fluid_mass_new = fluid_mass_node

    index_node_old = index_node

    for i in range(0, node_number):
        index_node_new[i] = index_node_old[i] + np.floor(fluid_evaporated_vol[i]/fluid_vol_unit)

    index_node_incre = index_node_new - index_node_old   
    index_node_cum = np.cumsum(index_node_incre) + index_node_old
    index_end = index_node_cum[-1]

    # index_node_incre[x+1 for x in range(0, node_number) if x == 0 ]
    index_node_incre_rev = index_node_incre + np.ones(node_number)

    for i in range(0, node_number):
        if index_node_incre[i] > 0:
            fluid_evaporated_vol_new[i] = 0

    results = [index_node_cum, fluid_mass_new, index_node_incre_rev, index_end, fluid_evaporated_vol_new]

    return results 

# pipe_node_new: index_node_cum
def Convert_vap_expansion(pipe_node_new, fluid_mass_node, node_delta, index, node_number, fluid_temp, fluid_specific_heat, fluid_density, fluid_quality, fluid_evaporated_vol, lng_density, ng_density, lng_specific_heat, ng_specific_heat, ambient_temp):

    index_here = node_number 
    pipe_properties = np.zeros((6, int(index_here)))
    pipe_properties[0,:] = ambient_temp

    # fluid_quality_tmp = 1- (1-fluid_quality[0])/node_delta[0]
    fluid_quality_tmp = fluid_quality[0]
    if fluid_quality_tmp > 1:
        fluid_quality_tmp = 1

    fluid_density_tmp = fluid_density[0]/node_delta[0]
    fluid_specific_heat_tmp = ng_specific_heat*fluid_quality_tmp + lng_specific_heat*(1- fluid_quality_tmp) 
    fluid_evaporated_tmp =  fluid_evaporated_vol[0]/node_delta[0]
    fluid_mass_tmp = fluid_mass_node[0]/node_delta[0]

    pipe_properties[:, :int(pipe_node_new[0])] = \
    np.transpose([[fluid_temp[0], fluid_specific_heat_tmp , fluid_density_tmp , fluid_quality_tmp, fluid_evaporated_vol[0], fluid_mass_tmp], ])
    flag_node[0] = 1

    for i in range(0, int(index_here)-1):
        if flag_node[i] == 0:
            pass
        else:
            index_cell_bef = int(pipe_node_new[i])
            if pipe_node_new[i] >= index_here:
                break
            index_cell_lat = int(pipe_node_new[i+1])
            if pipe_node_new[i+1] >= index_here:
                pipe_node_new[i+1] = index_here

        # fluid_quality_tmp = 1- (1-fluid_quality[i])/node_delta[i]
            fluid_quality_tmp = fluid_quality[i]
            if fluid_quality_tmp > 1:
                fluid_quality_tmp = 1
            if node_delta[i] == 0:
                node_delta[i] = 1
            fluid_density_tmp = fluid_density[i]/node_delta[i]
            fluid_specific_heat_tmp = ng_specific_heat*fluid_quality_tmp + lng_specific_heat*(1- fluid_quality_tmp) 
            fluid_evaporated_tmp =  fluid_evaporated_vol[i]/node_delta[i]
            fluid_mass_tmp = fluid_mass_node[i]/node_delta[i]           

            pipe_properties[:, index_cell_bef:index_cell_lat] =  \
            np.transpose([[fluid_temp[i], fluid_specific_heat_tmp, fluid_density_tmp, fluid_quality_tmp,  fluid_evaporated_tmp, fluid_mass_tmp],])          
            flag_node[index_cell_bef:index_cell_lat] = 1

            if pipe_node_new[i+1] >= index_here:
                break

    results = pipe_properties

    return results

# In[37]:

def Update_flow(fluid_temp_node, fluid_mass_node, fluid_specific_heat_node,fluid_density_node,fluid_quality_node, fluid_evaporated_vol, flow_rate, fluid_temp_ini, specific_heat_ini, density_ini, vol_unit, flag_node):
### Upgrading fluid flow

    fluid_temp_new = np.zeros_like(fluid_temp_node)
    fluid_specific_heat_new = np.zeros_like(fluid_temp_node)    
    fluid_density_new = np.zeros_like(fluid_temp_node)
    fluid_quality_new = np.zeros_like(fluid_temp_node)
    fluid_quality_new = np.zeros_like(fluid_temp_node)
    fluid_evaporated_vol_new = np.zeros_like(fluid_temp_node)
    fluid_mass_new = fluid_mass_node

    if flow_rate > 0:
        fluid_temp_new[1:] = fluid_temp_node[:-1]
        fluid_specific_heat_new[1:] = fluid_specific_heat_node[:-1]
        fluid_density_new[1:] = fluid_density_node[:-1]
        fluid_quality_new[1:] = fluid_quality_node[:-1]
        fluid_evaporated_vol_new[1:] = fluid_evaporated_vol[:-1]
        fluid_mass_new[1:] = fluid_mass_node[:-1]

        fluid_temp_new[0] = fluid_temp_ini
        fluid_specific_heat_new[0] = specific_heat_ini
        fluid_density_new[0] = density_ini
        fluid_quality_new[0] = 0
        fluid_evaporated_vol_new[0] = 0
        fluid_mass_new[0] = density_ini * vol_unit

    else:
        fluid_temp_new = fluid_temp_node
        fluid_specific_heat_new = fluid_specific_heat_node
        fluid_density_new = fluid_density_node
        fluid_quality_new = fluid_quality_node
        fluid_evaporated_vol_new = fluid_evaporated_vol
        fluid_mass_new = fluid_mass_node

    flag_node[1:] = flag_node[:-1]
    flag_node[0] = 1

    results = [fluid_temp_new, fluid_mass_new, fluid_specific_heat_new, fluid_density_new, fluid_quality_new, fluid_evaporated_vol_new, flag_node]

    return results

# In[39]:
### Time Integration
for T in range(0, time_number-1):

    if T == 21:      
        aa = 3
    heat_flow_ex[:,T], heatflow_ex[T] = Heat_trasnfer_ex(pipe_temp_node[:,T], ambient_temp_node, hc_ex, area_unit_ex, delta_time)
    heat_flow_in[:,T], heatflow_in[T] = Heat_trasnfer_in(fluid_temp_node[:,T], pipe_temp_node[:,T], fluid_quality_node[:,T], hc_in, area_unit_in, delta_time, node_number)

    pipe_specific_heat_node[:,T], pipe_heat_capa_node[:,T] = Update_pipe_pro(pipe_density, pipe_temp_node[:,T], pipe_volume_unit, node_number)
    pipe_temp_node[:,T+1] = Upgrade_pipe_temp(heat_flow_ex[:,T], heat_flow_in[:,T], pipe_heat_capa_node[:,T], pipe_temp_node[:,T], node_number)

    fluid_temp_node[:,T+1], fluid_mass_node[:,T], fluid_specific_heat_node[:,T+1], fluid_density_node[:,T+1], fluid_quality_node[:,T+1], fluid_evaporated_vol[:,T+1] =  \
    Upgrade_fluid_properties(heat_flow_in[:,T], fluid_mass_node[:,T], fluid_temp_node[:,T], fluid_quality_node[:,T], fluid_density_node[:,T], fluid_specific_heat_node[:,T], fluid_evaporated_vol[:,T], pressure, fluid_vol_unit, flag_node, lng_density, ng_density, lng_specific_heat, ng_specific_heat)

    pipe_node_new[:,T], fluid_mass_node[:,T], index_node_incre[:,T], index_end_variable, fluid_evaporated_vol[:,T+1] = \
    Vap_expansion(fluid_quality_node[:,T+1], fluid_mass_node[:,T], fluid_evaporated_vol[:,T+1], node_number, ng_density, lng_density, fluid_vol_unit, index_node)

    pipe_properties = Convert_vap_expansion(pipe_node_new[:,T], fluid_mass_node[:,T], index_node_incre[:,T], index_end_variable, node_number, fluid_temp_node[:,T+1], fluid_specific_heat_node[:,T+1], \
    fluid_density_node[:,T+1], fluid_quality_node[:,T+1], fluid_evaporated_vol[:,T+1], lng_density, ng_density, lng_specific_heat, ng_specific_heat, ambient_temperature)

    fluid_temp_node[:,T+1] = pipe_properties[0,:]
    fluid_specific_heat_node[:,T+1] = pipe_properties[1,:]
    fluid_density_node[:,T+1] = pipe_properties[2,:]
    fluid_quality_node[:,T+1] = pipe_properties[3,:]
    fluid_evaporated_vol[:,T+1] = pipe_properties[4,:]
    fluid_mass_node[:,T+1] = pipe_properties[5,:]

    fluid_temp_node[:,T+1],fluid_mass_node[:,T+1], fluid_specific_heat_node[:,T+1],fluid_density_node[:,T+1],fluid_quality_node[:,T+1], fluid_evaporated_vol[:,T+1], flag_node = \
    Update_flow(fluid_temp_node[:,T+1], fluid_mass_node[:,T+1], fluid_specific_heat_node[:,T+1],fluid_density_node[:,T+1],fluid_quality_node[:,T+1], fluid_evaporated_vol[:,T+1], \
    flow_rate[T], fluid_temp_ini, specific_heat_ini, density_ini,  fluid_vol_unit, flag_node)

    fluid_mass_node2[:,T+1] = fluid_density_node[:,T+1]*fluid_vol_unit

    index_history[T] = index_end_variable 
    masstot[T] = sum(fluid_density_node[:,T])*fluid_vol_unit
    # enthaltot[T] = sum(fluid_enthalpy_node[:,T])*fluid_volume_unit
    print ('')
    print(f'mass total: {masstot[T]:.2f}')
    print(f'enthal total: {enthaltot[T]:.2f}')
    print (' Iterate: %d     Time:        %.1f min.    Index history: %d ' %(T, times[T]/60, index_history[T] ))
    print ('')

# In[ ]:
N = int(np.floor(T/4))    

fig1 = plt.figure(figsize=(10,2.5), dpi=100)
plt.plot(distance_node[1:], pipe_temp_node[1:,N]-273.15, '-r', label= 'Pipe temp.1 ')
plt.plot(distance_node[1:], pipe_temp_node[1:,2*N]-273.15, '-b', label= 'Pipe temp.2')
plt.plot(distance_node[1:], pipe_temp_node[1:,3*N]-273.15, '-k', label= 'Pipe temp.3')
plt.plot(distance_node[1:], pipe_temp_node[1:,4*N]-273.15, '-y', label= 'Pipe temp.4')
plt.xlabel("Distance [m]")
plt.ylabel("Temperature [C]")
plt.grid(True)
plt.legend(loc='upper right')
plt.title('Pipe Temperature [C]')

fig2 = plt.figure(figsize=(10,2.5), dpi=100)
plt.plot(distance_node[1:], fluid_temp_node[1:,N]-273.15, '-r', label= 'Fluid temp.1')
plt.plot(distance_node[1:], fluid_temp_node[1:,2*N]-273.15, '-b', label= 'Fluid temp.2')
plt.plot(distance_node[1:], fluid_temp_node[1:,3*N]-273.15, '-k', label= 'Fluid temp.3')
plt.plot(distance_node[1:], fluid_temp_node[1:,4*N]-273.15, '-y', label= 'Fluid temp.4')
plt.xlabel("Distance [m]")
plt.ylabel("Temperature [C]")
plt.grid(True)
plt.legend(loc='upper right')
plt.title('fluid Temperature [C]')

fig22 = plt.figure(figsize=(10,2.5), dpi=100)
plt.plot(distance_node[1:], fluid_temp_node[1:,T]-273.15, '-r', label= 'Fluid temp.')
plt.xlabel("Distance [m]")
plt.ylabel("Temperature [C]")
plt.grid(True)
plt.legend(loc='upper right')
plt.title('fluid Temperature [C]')

fig22 = plt.figure(figsize=(10,2.5), dpi=100)
plt.plot(times[:T], masstot[:T], '-r', label= 'Fluid temp.')
plt.xlabel("Time [s]")
plt.ylabel("mass tot [kg]")
plt.grid(True)
plt.title('check mass conservation')

fig3 = plt.figure(figsize=(10,2.5), dpi=100)
plt.plot(distance_node[1:], fluid_quality_node[1:,N], label= 'Fluid quality 1')
plt.plot(distance_node[1:], fluid_quality_node[1:,2*N], label= 'Fluid quality 2')
plt.plot(distance_node[1:], fluid_quality_node[1:,3*N], label= 'Fluid quality 3')
plt.xlabel("Distance [m]")
plt.ylabel("Quality")
plt.grid(True)
plt.legend(loc='upper right')
plt.title('fluid Quality')

plt.title('pipe Temperature [K]')

### !!!!!!!!!  ###
fig57 = plt.figure(figsize=(10,5), dpi=100)
plt.plot(times[1:T]/60, heatflow_in[1:T]/3600, label= 'Pipe temperature')
plt.xlabel("Times [min]")
plt.ylabel("heat flux into LNG [W]")
plt.grid(True)
plt.legend(loc='upper right')
plt.title('Heat flow')
plt.xlim(0, 2)
### !!!!!!!!!  ###
# ss

# for i in range(0, T, 10):
#     fig6 = plt.figure(figsize=(10,6), dpi=100) 
#     # np.delete(H_nodes[:,T],IndexLNG+1) 
#     # np.delete(TempFluid_act2[:,i], IndexLNG+1) 
#     plt.plot(distance_node[:],  fluid_temp_node[:,i]-273.15, 'ro') 
#     plt.xlabel("Distance [m]")
#     plt.ylabel("Temperature [C]")
#     plt.grid(True)
#     plt.title('Temperature of pipe & fluid')
#     plt.ylim(-160, 35)

# for i in range(0, T, 1):
#     fig6 = plt.figure(figsize=(10,6), dpi=100) 
#     plt.plot(distance_node[:],  fluid_mass_node[:,i], 'ro') 
#     plt.xlabel("Distance [m]")
#     plt.ylabel("mass [kg]")
#     plt.grid(True)
#     plt.ylim(0, 0.1)

# ts = 2      # ts: time scale
# font1 = {'family': 'Times New Roman',
#          'color': 'blue',
#          'weight': 'normal',
#          'size': 20}
# images_final = []

# for i in range(0,120):
#     fig6 = plt.figure(figsize=(10,6), dpi=100)
#     plt.plot(distance_node[1:-1], pipe_temp_node[1:-1,ts*i]-273.15, label= 'Pipe temperature')
#     plt.plot(distance_node[1:-1], fluid_temp_node[1:-1,ts*i]-273.15, label= 'Fluid temperature')
#     plt.xlabel("Distance [m]")
#     plt.ylabel("Temperature [C]")
#     plt.grid(True)
#     plt.legend(loc='upper right')
#     plt.title('Temperature of pipe & fluid')
#     plt.ylim(-160, 35)
#     plt.text(0.2, 15, 'Time [min] :  ' + str(round(times[ts*i]/60,1)), fontdict= font1)

#     filename = openfiledir + '\\' + 'Fig' + '\\' + str(i)
#     plt.savefig(filename)

#     filename2 = filename + '.png'

#     images = Image.open(filename2) 
#     images_final.append(images)     

# Im = images_final[0]

# filename3 = filename + '0.30.gif'
# Im.save(filename3, save_all=True, append_images=images_final[1:], duration=100, loop=0) 

# for i in range(0,40):
#     fig6, ax1 = plt.subplots(figsize=(10,6), dpi=100)
#     ax1.plot(distance_node[1:], fluid_quality_node[1:,ts*i], label= 'fluid quality')
#     ax1.set_xlabel("Distance [m]")
#     ax1.set_ylabel("Quality")
#     ax2 = ax1.twinx() # instantiate a second axes that shares the same x-axis color_2 = 'tab:red' ax2.set_ylabel('Y2 value (red)', fontsize=14, color=color_2) ax2.plot(df.index, df.y2, marker='o', color=color_2)
#     ax2.plot(distance_node[1:], heat_flow_in[1:,ts*i], 'r-', label= 'heat flow')
#     ax1.set_ylim(0,1)
#     ax2.set_ylim(0,60)
#     ax2.set_ylabel("Heat Flow")
#     plt.grid(True)
#     plt.text(0.2, 15, 'Time [min] :  ' + str(round(times[ts*i]/60,1)), fontdict= font1)
JungHulk commented 2 years ago 
#!/usr/bin/env python
# coding: utf-8

# # Calculation of Pipe Heat transfer
# ## Consider pipe cooling and LNG evaporation
# 

# In[1]:

import os; os.environ['RPPREFIX'] = 'C:\\Program Files (x86)\\REFPROP'
from ctREFPROP.ctREFPROP import REFPROPFunctionLibrary
RP = REFPROPFunctionLibrary(os.environ['RPPREFIX'])
RP.SETPATHdll(os.environ['RPPREFIX'])

# In[2]:

from win32com.client import Dispatch
import os
import time; xxstart = time.time()

# In[3]:

import math
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image

openfiledir = os.getcwd()  # 현재 폴더 위치 

# In[4]:

MASS_BASE_SI = RP.GETENUMdll(0, "MASS BASE SI").iEnum   # output: 21

# In[8]:

lng_density = 430                     # [kg/m3]
lng_latent_heat = 510000              # [J/kg]

time_release = 1.5                      # [s]
time_total = 20                      # [s]
delta_time = 0.1                      # [s]

time_number = int(time_total/delta_time)         # [#].
times = np.linspace(0,time_total, time_number)

inner_diameter = 0.1143                 # [m]  50A
pipe_thickness = 0.003048                 # [m]  Sch. 10
# inner_diameter = 0.24448                 # [m]  50A
# pipe_thickness = 0.008687                 # [m]  Sch. 10
outer_diameter = inner_diameter + 2*pipe_thickness                                          # [m]
cross_area = 0.25*math.pi*inner_diameter**2                                                 # [m2]

# 2.53 kg, 36 초간 release
flow_rate = np.zeros(time_number) 
masstot = np.zeros(time_number) 
enthaltot = np.zeros(time_number) 
# flow_rate[:int(time_release/delta_time)] = 4535/lng_density/time_release/3600               # [m3/s]
flow_rate[:int(time_release/delta_time)] = 34790/lng_density/3600               # [m3/s]
# flow_rate[:int(time_total/delta_time)] = 2.53/lng_density/time_release                    # [m3/s]
flow_velocity = flow_rate[0]/cross_area
distance_progress = flow_velocity * delta_time

# In[9]:

delta_length = distance_progress                   # [m]
node_number = 250    #2000

line_length =  delta_length * node_number          # [m]
distance_node = np.linspace(0, line_length, node_number)

index_end_variable = 0

# In[10]:

ambient_temperature = 30 + 273.15                  # [K]
ambient_temp_node = ambient_temperature*np.ones(node_number)

# In[11]:

# In[12]:

pipe_volume_unit = 0.25*math.pi*(outer_diameter**2 - inner_diameter**2)*delta_length     # [m3]

inner_surface_tot = math.pi*inner_diameter*line_length                                   # [m2]
area_unit_in = math.pi*inner_diameter*delta_length                                       # [m2]

fluid_vol_unit = cross_area * delta_length
inner_surface_tot, area_unit_in, cross_area

# In[13]:

# In[14]:

# pipe material: SUS316L
pipe_thermal_conductivity = 16                     # [W/mk]
pipe_density = 8030                                # [kg/m3]
pipe_specific_heat = 3.1121*ambient_temperature - 0.0051*ambient_temperature**2                          # [J/kg K]

# In[15]:

pipe_heat_capa_unit = pipe_specific_heat*pipe_density*pipe_volume_unit                   # [J/K]
pipe_heat_capa_unit

# In[16]:

insulation_thickness = 0.03                                                                               # [m]
insulation_outer_diameter = outer_diameter + 2 * pipe_thickness                                           # [m]

insulation_volume_unit = 0.25*math.pi*(insulation_outer_diameter**2 - outer_diameter**2)*delta_length     # [m3]
insulation_volume_unit, pipe_thickness

area_unit_ex = math.pi*insulation_outer_diameter*delta_length 

# In[17]:

# Insulation material: Rock Wool
insulation_density = 55                                                                                   # [kg/m3]
insulation_specific_heat = 840                                                                            # [J/kg K]
insulation_thermal_conductivity = 0.03                                                                    # [W/mK]

insulation_heat_capa_unit = insulation_specific_heat*insulation_density*insulation_volume_unit            # [J/K]
insulation_heat_capa_unit

# In[18]:

hc_in= 800                                                                      # [W/m2K]
hc_ex = 8                                                                      # [W/m2K]
pipe_heat_coeff = pipe_thermal_conductivity/pipe_thickness                                # [W/m2K]
insulation_heat_coeff = insulation_thermal_conductivity/insulation_thickness              # [W/m2K]

# overall_heat_coeff = 1/(1/inner_heat_coeff + 1/pipe_heat_coeff + 1/insulation_heat_coeff +1/outer_heat_coeff)
overall_heat_coeff = 1/(1/hc_in + 1/pipe_heat_coeff +1/hc_ex)
overall_heat_coeff

# In[19]:

################### fluid properties #########################
fluid_temp_ini = -120 + 273.15                                                      # [K]
pressure = (16.2 + 1)* 101300                                                               # [Pa]

# In[22]:

fluid = "methane"
znit = [1.]
flag_node = np.zeros(node_number)

lng_density = RP.REFPROPdll(fluid, "PT", "D", MASS_BASE_SI, 1, 0, pressure, fluid_temp_ini, znit).Output[0]         # [kg/m3]
lng_specific_heat = RP.REFPROPdll(fluid, "PT", "Cp", MASS_BASE_SI, 1, 0, pressure, fluid_temp_ini, znit).Output[0]  # [J/kg K]

ng_density = RP.REFPROPdll(fluid, "PQ", "D", MASS_BASE_SI, 1, 0, pressure, 1, znit).Output[0]         # [kg/m3]
ng_specific_heat = RP.REFPROPdll(fluid, "PQ", "Cp", MASS_BASE_SI, 1, 0, pressure, 1, znit).Output[0]  # [J/kg K]

specific_heat_ini = lng_specific_heat
density_ini = lng_density 

ng_density, ng_specific_heat, lng_density, lng_specific_heat
index_node = np.linspace(0,node_number-1,node_number)
index_node_incre= np.zeros((node_number, time_number))

# In[23]:

index_history = np.zeros(time_number)
heatflow_ex = np.zeros(time_number)
heatflow_in = np.zeros(time_number)

passing_vapor = np.zeros(time_number)
passing_liquid = np.zeros(time_number)

pipe_temp_node = np.zeros((node_number, time_number))
pipe_specific_heat_node = np.zeros((node_number, time_number))
pipe_heat_capa_node = np.zeros((node_number, time_number))

insulation_temp_node = np.zeros((node_number, time_number))
fluid_temp_node = np.zeros((node_number, time_number))
fluid_quality_node = np.zeros((node_number, time_number))
fluid_density_node = np.zeros((node_number, time_number))
fluid_mass_node = np.zeros((node_number, time_number))
fluid_mass_node2 = np.zeros((node_number, time_number))
fluid_specific_heat_node = np.zeros((node_number, time_number))
fluid_evaporated_vol = np.zeros((node_number, time_number))

heat_flow_ex = np.zeros((node_number, time_number))
heat_flow_in = np.zeros((node_number, time_number))

heat_flux = np.zeros((node_number, time_number))
pipe_node_new = np.zeros((node_number, time_number))

## Initialization
pipe_temp_node[:,0] = ambient_temperature                                                    # [K]
fluid_temp_node[:,:] = ambient_temperature  
fluid_quality_node[:,0] = 0  
fluid_density_node[:,0] = 0
fluid_specific_heat_node[:,0] = lng_specific_heat

# insulation_temp_node[:] = ambient_temperature                                              # [K]

# In[24]:

def Heat_trasnfer_ex(pipe_temp_node, ambient_temp_node, hc_ex, area_unit_ex, delta_time):
    heat_flow_ex = hc_ex*area_unit_ex*(ambient_temp_node - pipe_temp_node) * delta_time

    heat_flow_ex_sum = np.sum(heat_flow_ex)

    results = [heat_flow_ex, heat_flow_ex_sum]
    return results

# In[25]:

def Heat_trasnfer_in(fluid_temp_node, pipe_temp_node, fluid_quality, hc_in, area_unit_in, delta_time, node_number):

    hc_in_node = np.zeros(node_number)
    hc_in_node[:] = hc_in

    for i in range(0, node_number):

        if pipe_temp_node[i] - fluid_temp_node[i] <5:
            hc_in_node[i] = 1500

        if fluid_quality[i] == 1:
            hc_in_node[i] = 20

    heat_flow_in = hc_in_node*area_unit_in*(pipe_temp_node - fluid_temp_node) * delta_time
    heat_flow_in_sum = np.sum(heat_flow_in)
    results = [heat_flow_in, heat_flow_in_sum]
    return results

def Update_pipe_pro(pipe_density, pipe_temp_node, pipe_volume_unit, node_number):

    pipe_specific_heat_new = np.zeros(node_number)
    pipe_heat_capa_new = np.zeros(node_number)

    pipe_specific_heat_new = 3.1121*pipe_temp_node - 0.0051*pipe_temp_node**2
    pipe_heat_capa_new = pipe_specific_heat_new*pipe_density*pipe_volume_unit

    results = [pipe_specific_heat_new, pipe_heat_capa_new]

    return results

# In[31]:

def Upgrade_pipe_temp(heat_flow_ex, heat_flow_in, pipe_specific_heat_node, pipe_temp_node, node_number):
    temp_delt = np.zeros(node_number)

    temp_delt = (heat_flow_ex - heat_flow_in)/pipe_specific_heat_node
    pipe_temp_node = pipe_temp_node + temp_delt

    return pipe_temp_node

# In[36]:

def Upgrade_fluid_properties(heat_flow_in, fluid_mass_node, fluid_temp_node, fluid_quality_node, fluid_density_node, fluid_specific_heat_node, fluid_evaporated_vol, pressure, fluid_vol_unit, flag_node, lng_density, ng_density, lng_specific_heat, ng_specific_heat):
# calculating the mass based fluid properties
    fluid_heat_capa = np.zeros(node_number)
    fluid_temp_new = fluid_temp_node
    fluid_specific_heat_new = np.zeros(node_number)
    fluid_density_new = np.zeros(node_number)
    fluid_mass_new  = np.zeros(node_number)
    fluid_quality_new = np.zeros(node_number)
    fluid_evaporated_vol_new = np.zeros(node_number)

    for i in range(0, node_number):
        mass_evap = 0
        # if i == 175:
        #     aa = 0
        if flag_node[i] >0:

            if fluid_quality_node[i] == 0:      # Liquid case

                mass_evap = heat_flow_in[i]/lng_latent_heat
                mass_tot = fluid_mass_node[i]

                fluid_evaporated_vol_new[i] = fluid_evaporated_vol[i] + mass_evap/ng_density        # [m3]

                if mass_tot ==0:
                    fluid_quality_new[i] = 0
                else:
                    fluid_quality_new[i] = mass_evap/mass_tot
                fluid_temp_new[i] = fluid_temp_node[i]

                fluid_specific_heat_new[i] = ng_specific_heat*fluid_quality_new[i] + lng_specific_heat*(1-fluid_quality_new[i])
                fluid_density_new[i] = fluid_mass_node[i]/fluid_vol_unit
                fluid_mass_new[i] = fluid_mass_node[i]

            elif fluid_quality_node[i] >= 1:      # Vapor case

                fluid_quality_node[i] = 1         
                # fluid_heat_capa[i] = fluid_density_node[i]*fluid_specific_heat_node[i]*fluid_vol_unit
                fluid_heat_capa[i] = ng_density*ng_specific_heat*fluid_vol_unit                

                if fluid_density_node[i] == 0:
                    fluid_density_node[i] = fluid_density_node[i-1]
                    fluid_heat_capa[i] = fluid_density_node[i]*fluid_specific_heat_node[i]*fluid_vol_unit                    
                if fluid_heat_capa[i] == 0:
                    print('error in line 293')

                delta_temp = heat_flow_in[i]/fluid_heat_capa[i]
                fluid_evaporated_vol_new[i] = fluid_evaporated_vol[i]

                fluid_temp_new[i] = fluid_temp_node[i] + delta_temp
                fluid_specific_heat_new[i] = ng_specific_heat
                fluid_density_new[i] = fluid_mass_node[i]/fluid_vol_unit
                fluid_quality_new[i] = fluid_quality_node[i]
                fluid_mass_new[i] = fluid_mass_node[i]

            elif fluid_quality_node[i] < 1 and  fluid_quality_node[i] > 0 :                               # Two phase case

                mass_evap = heat_flow_in[i]/lng_latent_heat
                mass_tot = fluid_mass_node[i]   

                fluid_evaporated_vol_new[i] = fluid_evaporated_vol[i] + mass_evap/ng_density
                if mass_tot ==0:
                    fluid_quality_delta = 0
                else:
                    fluid_quality_delta = mass_evap/mass_tot
                fluid_quality_new[i] = fluid_quality_node[i] + fluid_quality_delta

                if fluid_quality_new[i] > 1:
                    fluid_quality_new[i] = 1

                fluid_density_new[i] = fluid_mass_node[i]/fluid_vol_unit 
                # if fluid_density_new[i] <= 6:
                #     fluid_density_new[i] = 6
                #     fluid_quality_new[i] = 1

                fluid_temp_new[i] = fluid_temp_node[i]            
                fluid_specific_heat_new[i] = ng_specific_heat*fluid_quality_new[i] + lng_specific_heat*(1-fluid_quality_new[i])     
                fluid_mass_new[i] = fluid_mass_node[i]
        else:
            fluid_temp_new[i] = fluid_temp_node[i]
            fluid_specific_heat_new[i] = fluid_specific_heat_node[i]
            fluid_density_new[i] = fluid_density_node[i]
            fluid_quality_new[i] = fluid_quality_node[i]
            fluid_evaporated_vol_new[i] = fluid_evaporated_vol[i]
            fluid_mass_new[i] = fluid_mass_node[i]         

    results = [fluid_temp_new, fluid_mass_new, fluid_specific_heat_new, fluid_density_new, fluid_quality_new, fluid_evaporated_vol_new] 

    return results

def Vap_expansion(fluid_quality_node, fluid_mass_node, fluid_evaporated_vol, node_number, ng_density, lng_density, fluid_vol_unit,  index_node):
## pipe_node_new: index_node_cum
    index_node_new = np.zeros(node_number)
    fluid_evaporated_vol_new = fluid_evaporated_vol
    index_node_incre = np.zeros(node_number)  
    fluid_mass_new = fluid_mass_node

    index_node_old = index_node

    for i in range(0, node_number):
        index_node_new[i] = index_node_old[i] + np.floor(fluid_evaporated_vol[i]/fluid_vol_unit)

    index_node_incre = index_node_new - index_node_old   
    index_node_cum = np.cumsum(index_node_incre) + index_node_old
    index_end = index_node_cum[-1]

    # index_node_incre[x+1 for x in range(0, node_number) if x == 0 ]
    index_node_incre_rev = index_node_incre + np.ones(node_number)

    for i in range(0, node_number):
        if index_node_incre[i] > 0:
            fluid_evaporated_vol_new[i] = fluid_evaporated_vol[i] - fluid_vol_unit

    results = [index_node_cum, fluid_mass_new, index_node_incre_rev, index_end, fluid_evaporated_vol_new]

    return results 

# pipe_node_new: index_node_cum
def Convert_vap_expansion(pipe_node, fluid_mass_node, node_delta, index, node_number, fluid_temp, fluid_specific_heat, fluid_density, fluid_quality, \
                          fluid_evaporated_vol, lng_density, ng_density, lng_specific_heat, ng_specific_heat, vol_unit, ambient_temp):

    index_here = node_number 
    pipe_properties = np.zeros((6, int(index_here)))
    pipe_properties[0,:] = ambient_temp

    # fluid_quality_tmp = 1- (1-fluid_quality[0])/node_delta[0]
    fluid_quality_tmp = fluid_quality[0]
    if fluid_quality_tmp > 1:
        fluid_quality_tmp = 1

    fluid_density_tmp = fluid_mass_node[0]/vol_unit/node_delta[0]
    fluid_specific_heat_tmp = ng_specific_heat*fluid_quality_tmp + lng_specific_heat*(1- fluid_quality_tmp) 
    fluid_evaporated_tmp =  fluid_evaporated_vol[0]/node_delta[0]
    fluid_mass_tmp = fluid_mass_node[0]/node_delta[0]

    pipe_properties[:, :int(pipe_node[0])+1] = \
    np.transpose([[fluid_temp[0], fluid_specific_heat_tmp , fluid_density_tmp , fluid_quality_tmp, fluid_evaporated_vol[0], fluid_mass_tmp], ])
    flag_node[0] = 1

    for i in range(int(index_here)-1):
        if flag_node[i] == 0:
            pass
        else:
            index_cell_bef = int(pipe_node[i]+1)
            if pipe_node[i] >= index_here:
                break
            index_cell_lat = int(pipe_node[i+1]+1)
            if index_cell_lat >= index_here:
                index_cell_lat = index_here

        # fluid_quality_tmp = 1- (1-fluid_quality[i])/node_delta[i]
            fluid_quality_tmp = fluid_quality[i]
            if fluid_quality_tmp >= 1:
                fluid_quality_tmp = 1
            if node_delta[i] == 0:
                node_delta[i] = 1
            fluid_density_tmp = fluid_mass_node[i+1]/vol_unit/node_delta[i+1]
            fluid_specific_heat_tmp = ng_specific_heat*fluid_quality_tmp + lng_specific_heat*(1- fluid_quality_tmp) 
            fluid_evaporated_tmp =  fluid_evaporated_vol[i+1]/node_delta[i+1]
            fluid_mass_tmp = fluid_mass_node[i+1]/node_delta[i+1]           

            pipe_properties[:, index_cell_bef:index_cell_lat] =  \
            np.transpose([[fluid_temp[i+1], fluid_specific_heat_tmp, fluid_density_tmp, fluid_quality_tmp,  fluid_evaporated_tmp, fluid_mass_tmp],])          
            flag_node[index_cell_bef:index_cell_lat] = 1

            if pipe_node[i+1] >= index_here:
                break

    results = pipe_properties

    return results

# In[37]:

def Update_flow(fluid_temp_node, fluid_mass_node, fluid_specific_heat_node,fluid_density_node,fluid_quality_node, fluid_evaporated_vol, flow_rate, fluid_temp_ini, specific_heat_ini, density_ini, vol_unit, flag_node):
### Upgrading fluid flow

    fluid_temp_new = np.zeros_like(fluid_temp_node)
    fluid_specific_heat_new = np.zeros_like(fluid_temp_node)    
    fluid_density_new = np.zeros_like(fluid_temp_node)
    fluid_quality_new = np.zeros_like(fluid_temp_node)
    fluid_quality_new = np.zeros_like(fluid_temp_node)
    fluid_evaporated_vol_new = np.zeros_like(fluid_temp_node)
    fluid_mass_new = fluid_mass_node

    if flow_rate > 0:
        fluid_temp_new[1:] = fluid_temp_node[:-1]
        fluid_specific_heat_new[1:] = fluid_specific_heat_node[:-1]
        fluid_density_new[1:] = fluid_density_node[:-1]
        fluid_quality_new[1:] = fluid_quality_node[:-1]
        fluid_evaporated_vol_new[1:] = fluid_evaporated_vol[:-1]
        fluid_mass_new[1:] = fluid_mass_node[:-1]

        fluid_temp_new[0] = fluid_temp_ini
        fluid_specific_heat_new[0] = specific_heat_ini
        fluid_density_new[0] = density_ini
        fluid_quality_new[0] = 0
        fluid_evaporated_vol_new[0] = 0
        fluid_mass_new[0] = density_ini * vol_unit

    else:
        fluid_temp_new = fluid_temp_node
        fluid_specific_heat_new = fluid_specific_heat_node
        fluid_density_new = fluid_density_node
        fluid_quality_new = fluid_quality_node
        fluid_evaporated_vol_new = fluid_evaporated_vol
        fluid_mass_new = fluid_mass_node

    flag_node[1:] = flag_node[:-1]
    flag_node[0] = 1

    results = [fluid_temp_new, fluid_mass_new, fluid_specific_heat_new, fluid_density_new, fluid_quality_new, fluid_evaporated_vol_new, flag_node]

    return results

# In[39]:
### Time Integration
for T in range(0, time_number-1):

    if T == 191:      
        aa = 3
    heat_flow_ex[:,T], heatflow_ex[T] = Heat_trasnfer_ex(pipe_temp_node[:,T], ambient_temp_node, hc_ex, area_unit_ex, delta_time)
    heat_flow_in[:,T], heatflow_in[T] = Heat_trasnfer_in(fluid_temp_node[:,T], pipe_temp_node[:,T], fluid_quality_node[:,T], hc_in, area_unit_in, delta_time, node_number)

    pipe_specific_heat_node[:,T], pipe_heat_capa_node[:,T] = Update_pipe_pro(pipe_density, pipe_temp_node[:,T], pipe_volume_unit, node_number)
    pipe_temp_node[:,T+1] = Upgrade_pipe_temp(heat_flow_ex[:,T], heat_flow_in[:,T], pipe_heat_capa_node[:,T], pipe_temp_node[:,T], node_number)

    fluid_temp_node[:,T+1], fluid_mass_node[:,T], fluid_specific_heat_node[:,T+1], fluid_density_node[:,T+1], fluid_quality_node[:,T+1], fluid_evaporated_vol[:,T+1] =  \
    Upgrade_fluid_properties(heat_flow_in[:,T], fluid_mass_node[:,T], fluid_temp_node[:,T], fluid_quality_node[:,T], fluid_density_node[:,T], fluid_specific_heat_node[:,T], fluid_evaporated_vol[:,T], pressure, fluid_vol_unit, flag_node, lng_density, ng_density, lng_specific_heat, ng_specific_heat)

    pipe_node_new[:,T], fluid_mass_node[:,T], index_node_incre[:,T], index_end_variable, fluid_evaporated_vol[:,T+1] = \
    Vap_expansion(fluid_quality_node[:,T+1], fluid_mass_node[:,T], fluid_evaporated_vol[:,T+1], node_number, ng_density, lng_density, fluid_vol_unit, index_node)

    pipe_properties = Convert_vap_expansion(pipe_node_new[:,T], fluid_mass_node[:,T], index_node_incre[:,T], index_end_variable, node_number, fluid_temp_node[:,T+1], fluid_specific_heat_node[:,T+1], \
    fluid_density_node[:,T+1], fluid_quality_node[:,T+1], fluid_evaporated_vol[:,T+1], lng_density, ng_density, lng_specific_heat, ng_specific_heat, fluid_vol_unit, ambient_temperature)

    fluid_temp_node[:,T+1] = pipe_properties[0,:]
    fluid_specific_heat_node[:,T+1] = pipe_properties[1,:]
    fluid_density_node[:,T+1] = pipe_properties[2,:]
    fluid_quality_node[:,T+1] = pipe_properties[3,:]
    fluid_evaporated_vol[:,T+1] = pipe_properties[4,:]
    fluid_mass_node[:,T+1] = pipe_properties[5,:]

    fluid_temp_node[:,T+1],fluid_mass_node[:,T+1], fluid_specific_heat_node[:,T+1],fluid_density_node[:,T+1],fluid_quality_node[:,T+1], fluid_evaporated_vol[:,T+1], flag_node = \
    Update_flow(fluid_temp_node[:,T+1], fluid_mass_node[:,T+1], fluid_specific_heat_node[:,T+1],fluid_density_node[:,T+1],fluid_quality_node[:,T+1], fluid_evaporated_vol[:,T+1], \
    flow_rate[T], fluid_temp_ini, specific_heat_ini, density_ini, fluid_vol_unit, flag_node)

    fluid_mass_node2[:,T+1] = fluid_density_node[:,T+1]*fluid_vol_unit

    index_history[T] = index_end_variable 
    masstot[T] = sum(fluid_mass_node[:,T])
    # enthaltot[T] = sum(fluid_enthalpy_node[:,T])*fluid_volume_unit
    print ('')
    print(f'mass total: {masstot[T]:.2f}')
    print(f'heat flux to Liq.: {heatflow_in[T]:.2f}')
    print (' Iterate: %d     Time:        %.1f min.    Index history: %d ' %(T, times[T]/60, index_history[T] ))
    print ('')

# In[ ]:
N = int(np.floor(T/4))    

fig1 = plt.figure(figsize=(10,2.5), dpi=100)
plt.plot(distance_node[1:], pipe_temp_node[1:,N]-273.15, '-r', label= 'Time: '+str(np.round(times[N]/60,2))+' Min.')
plt.plot(distance_node[1:], pipe_temp_node[1:,2*N]-273.15, '-b', label= 'Time: '+str(np.round(times[2*N]/60,2))+' Min.')
plt.plot(distance_node[1:], pipe_temp_node[1:,3*N]-273.15, '-k', label= 'Time: '+str(np.round(times[3*N]/60,2))+' Min.')
plt.plot(distance_node[1:], pipe_temp_node[1:,4*N]-273.15, '-y', label= 'Time: '+str(np.round(times[4*N]/60,2))+' Min.')
plt.xlabel("Distance [m]")
plt.ylabel("Temperature [C]")
plt.grid(True)
plt.legend(loc='upper right')
plt.title('Pipe Temperature [C]')

fig2 = plt.figure(figsize=(10,2.5), dpi=100)
plt.plot(distance_node[1:], fluid_temp_node[1:,N]-273.15, '-r', label= 'Time: '+str(np.round(times[N]/60,2))+' Min.')
plt.plot(distance_node[1:], fluid_temp_node[1:,2*N]-273.15, '-b', label= 'Time: '+str(np.round(times[2*N]/60,2))+' Min.')
plt.plot(distance_node[1:], fluid_temp_node[1:,3*N]-273.15, '-k', label= 'Time: '+str(np.round(times[3*N]/60,2))+' Min.')
plt.plot(distance_node[1:], fluid_temp_node[1:,4*N]-273.15, '-y', label= 'Time: '+str(np.round(times[4*N]/60,2))+' Min.')
plt.xlabel("Distance [m]")
plt.ylabel("Temperature [C]")
plt.grid(True)
plt.legend(loc='upper right')
plt.title('fluid Temperature [C]')

fig22 = plt.figure(figsize=(10,2.5), dpi=100)
plt.plot(distance_node[1:], fluid_temp_node[1:,T]-273.15, '-r', label= 'Fluid temp.')
plt.xlabel("Distance [m]")
plt.ylabel("Temperature [C]")
plt.grid(True)
plt.legend(loc='upper right')
plt.title('fluid Temperature [C]')

fig22 = plt.figure(figsize=(10,2.5), dpi=100)
plt.plot(times[:T], masstot[:T], '-r', label= 'Fluid temp.')
plt.xlabel("Time [s]")
plt.ylabel("mass tot [kg]")
plt.grid(True)
plt.title('check mass conservation')

fig3 = plt.figure(figsize=(10,2.5), dpi=100)
plt.plot(distance_node[1:], fluid_quality_node[1:,N], label= 'Time: '+str(np.round(times[N]/60,2))+' Min.')
plt.plot(distance_node[1:], fluid_quality_node[1:,2*N], label= 'Time: '+str(np.round(times[2*N]/60,2))+' Min.')
plt.plot(distance_node[1:], fluid_quality_node[1:,3*N], label= 'Time: '+str(np.round(times[3*N]/60,2))+' Min.')
plt.plot(distance_node[1:], fluid_quality_node[1:,4*N], label= 'Time: '+str(np.round(times[4*N]/60,2))+' Min.')
plt.xlabel("Distance [m]")
plt.ylabel("Quality")
plt.grid(True)
plt.legend(loc='upper right')
plt.title('fluid Quality')

### !!!!!!!!!  ###
fig57 = plt.figure(figsize=(10,5), dpi=100)
plt.plot(times[1:T]/60, heatflow_in[1:T]/3600, label= 'Heat flow into pipe')
plt.xlabel("Times [min]")
plt.ylabel("heat flux into LNG [W]")
plt.grid(True)
plt.legend(loc='upper right')
plt.title('Heat flow')
plt.xlim(0, 2)
### !!!!!!!!!  ###
# ss

# for i in range(0, T, 10):
#     fig6 = plt.figure(figsize=(10,6), dpi=100) 
#     # np.delete(H_nodes[:,T],IndexLNG+1) 
#     # np.delete(TempFluid_act2[:,i], IndexLNG+1) 
#     plt.plot(distance_node[:],  fluid_temp_node[:,i]-273.15, 'ro') 
#     plt.xlabel("Distance [m]")
#     plt.ylabel("Temperature [C]")
#     plt.grid(True)
#     plt.title('Temperature of pipe & fluid')
#     plt.ylim(-160, 35)

# for i in range(0, T, 1):
#     fig6 = plt.figure(figsize=(10,6), dpi=100) 
#     plt.plot(distance_node[:],  fluid_mass_node[:,i], 'ro') 
#     plt.xlabel("Distance [m]")
#     plt.ylabel("mass [kg]")
#     plt.grid(True)
#     plt.ylim(0, 0.1)

# ts = 2      # ts: time scale
# font1 = {'family': 'Times New Roman',
#          'color': 'blue',
#          'weight': 'normal',
#          'size': 20}
# images_final = []

# for i in range(0,120):
#     fig6 = plt.figure(figsize=(10,6), dpi=100)
#     plt.plot(distance_node[1:-1], pipe_temp_node[1:-1,ts*i]-273.15, label= 'Pipe temperature')
#     plt.plot(distance_node[1:-1], fluid_temp_node[1:-1,ts*i]-273.15, label= 'Fluid temperature')
#     plt.xlabel("Distance [m]")
#     plt.ylabel("Temperature [C]")
#     plt.grid(True)
#     plt.legend(loc='upper right')
#     plt.title('Temperature of pipe & fluid')
#     plt.ylim(-160, 35)
#     plt.text(0.2, 15, 'Time [min] :  ' + str(round(times[ts*i]/60,1)), fontdict= font1)

#     filename = openfiledir + '\\' + 'Fig' + '\\' + str(i)
#     plt.savefig(filename)

#     filename2 = filename + '.png'

#     images = Image.open(filename2) 
#     images_final.append(images)     

# Im = images_final[0]

# filename3 = filename + '0.30.gif'
# Im.save(filename3, save_all=True, append_images=images_final[1:], duration=100, loop=0) 

# for i in range(0,40):
#     fig6, ax1 = plt.subplots(figsize=(10,6), dpi=100)
#     ax1.plot(distance_node[1:], fluid_quality_node[1:,ts*i], label= 'fluid quality')
#     ax1.set_xlabel("Distance [m]")
#     ax1.set_ylabel("Quality")
#     ax2 = ax1.twinx() # instantiate a second axes that shares the same x-axis color_2 = 'tab:red' ax2.set_ylabel('Y2 value (red)', fontsize=14, color=color_2) ax2.plot(df.index, df.y2, marker='o', color=color_2)
#     ax2.plot(distance_node[1:], heat_flow_in[1:,ts*i], 'r-', label= 'heat flow')
#     ax1.set_ylim(0,1)
#     ax2.set_ylim(0,60)
#     ax2.set_ylabel("Heat Flow")
#     plt.grid(True)
#     plt.text(0.2, 15, 'Time [min] :  ' + str(round(times[ts*i]/60,1)), fontdict= font1)