Basic HTE

[JungHulk]

Tank_information

# -*- coding: utf-8 -*-
"""
Created on Wed Jul 20 09:37:03 2022

@author: 252914
"""
import numpy as np

class Cargo_tank():

    # def Tank_info(Tank_type): # 0: Membrane, 1: Flex, 2: Type-B, 3:Type-C

    Type_tank_list = ['Membrane','Membrane-Flex', 'Type-B', 'Type-C']
    Type_tank = Type_tank_list[0]
    if Type_tank == Type_tank_list[3]:
        Circle_ratio = 3
    else:
        Circle_ratio = 1.7

    Model = "SN2434"

    if Model == "7_5k_TypeC":

        BOR = 0.2           # %/day
        FR = 90             # %    
        ## Tank 정보는 Model 표 만들어서 거기서 가져오는것으로 만들기
        Vol_tank_1 = 3750         # m3
        Vol_tank_2 = 0        # m3
        Vol_tank_3 = 0    # m3
        Vol_tank_4 = 0    # m3
        Vol_tank_list = np.array([Vol_tank_1, Vol_tank_2, Vol_tank_3, Vol_tank_4])
        Vol_tank = np.sum(Vol_tank_list)

        Area_tank_1 = 2874/2         # m2
        Area_tank_2 = 0         # m2
        Area_tank_3 = 0  # m2
        Area_tank_4 = 0  # m2
        Area_tank_list = np.array([Area_tank_1, Area_tank_2, Area_tank_3, Area_tank_4])
        Area_tank = np.sum(Area_tank_list)

        Tank_info = [Type_tank, Circle_ratio, Vol_tank_list, Vol_tank, Area_tank_list, Area_tank, BOR, FR]

        Ins_thick = 300            # mm     
        Wall_thick = 24            # mm 
        M_wall = 550000/2            # kg SUS weight
        Cp_wall = 0.35             # kJ/kg/C
        Cond_wall = 17             # kJ/h/m/C

        M_ins = 44000/2            # kg Ins weight
        Cp_ins = 1.071             # kJ/kg/C
        D_ins = 118                # kg/m3
        Cond_ins = 0.0837           # kJ/h/m/C

        ht_NG = 37.46              # kJ/h/m2/C
        ht_atm = 19.5              # kJ/h/m2/C
        Tank_prop = [Ins_thick, M_wall, Cp_wall, M_ins, Cp_ins, D_ins, Cond_ins, ht_NG, ht_atm]

        sec = 0  # secondary barrier location in insulation

    elif Model == "SN2430":

        BOR = 0.239           # %/day
        FR = 90             # %    
        ## Tank 정보는 Model 표 만들어서 거기서 가져오는것으로 만들기
        Vol_tank_1 = 3777         # m3
        Vol_tank_2 = 3777        # m3
        Vol_tank_3 = 0    # m3
        Vol_tank_4 = 0    # m3
        Vol_tank_list = np.array([Vol_tank_1, Vol_tank_2, Vol_tank_3, Vol_tank_4])
        Vol_tank = np.sum(Vol_tank_list)

        Area_tank_1 = 2874 /2        # m2
        Area_tank_2 = 2874 /2        # m2
        Area_tank_3 = 0  # m2
        Area_tank_4 = 0  # m2
        Area_tank_list = np.array([Area_tank_1, Area_tank_2, Area_tank_3, Area_tank_4])
        Area_tank = np.sum(Area_tank_list)

        Tank_info = [Type_tank, Circle_ratio, Vol_tank_list, Vol_tank, Area_tank_list, Area_tank, BOR, FR]

        Ins_thick = 300            # mm     
        Wall_thick = 24            # mm 
        M_wall = 550000            # kg SUS weight
        Cp_wall = 0.35             # kJ/kg/C
        Cond_wall = 17             # kJ/h/m/C

        M_ins = 44000            # kg Ins weight
        Cp_ins = 1.071             # kJ/kg/C
        D_ins = 40                # kg/m3
        Cond_ins = 0.0837           # kJ/h/m/C

        ht_NG = 420              # kJ/h/m2/C
        ht_atm = 50              # kJ/h/m2/C
        Tank_prop = [Ins_thick, M_wall, Cp_wall, M_ins, Cp_ins, D_ins, Cond_ins, ht_NG, ht_atm]

        sec = 0  # secondary barrier location in insulation

    elif Model == "6k_Mem": # SN2586

        BOR = 0.345           # %/day
        FR = 98.5             # %    
        ## Tank 정보는 Model 표 만들어서 거기서 가져오는것으로 만들기
        Vol_tank_1 = 6000         # m3
        Vol_tank_2 = 0        # m3
        Vol_tank_3 = 0    # m3
        Vol_tank_4 = 0    # m3
        Vol_tank_list = np.array([Vol_tank_1, Vol_tank_2, Vol_tank_3, Vol_tank_4])
        Vol_tank = np.sum(Vol_tank_list)

        Area_tank_1 = 2075         # m2
        Area_tank_2 = 0         # m2
        Area_tank_3 = 0  # m2
        Area_tank_4 = 0  # m2
        Area_tank_list = np.array([Area_tank_1, Area_tank_2, Area_tank_3, Area_tank_4])
        Area_tank = np.sum(Area_tank_list)

        Tank_info = [Type_tank, Circle_ratio, Vol_tank_list, Vol_tank, Area_tank_list, Area_tank, BOR, FR]

        Ins_thick = 270            # mm    
        Wall_thick = 1.2            # mm 
        M_wall = 23912            # kg SUS weight
        Cp_wall = 0.46             # kJ/kg/C
        Cond_wall = 35             # kJ/h/m/C

        M_ins = 175343            # kg Ins weight
        Cp_ins = 1.071             # kJ/kg/C
        D_ins = 118                # kg/m3
        Cond_ins = 0.0837           # kJ/h/m/C

        ht_NG = 37.76              # kJ/h/m2/C
        ht_atm = 24              # kJ/h/m2/C
        Tank_prop = [Ins_thick, M_wall, Cp_wall, M_ins, Cp_ins, D_ins, Cond_ins, ht_NG, ht_atm]

        sec = 0.1  # secondary barrier location in insulation

    elif Model == "12k_Type-B": # SN2434

        BOR = 0.2           # %/day
        FR = 98.0            # %    
        ## Tank 정보는 Model 표 만들어서 거기서 가져오는것으로 만들기
        Vol_tank_1 = 12000         # m3
        Vol_tank_2 = 0        # m3
        Vol_tank_3 = 0    # m3
        Vol_tank_4 = 0    # m3
        Vol_tank_list = np.array([Vol_tank_1, Vol_tank_2, Vol_tank_3, Vol_tank_4])
        Vol_tank = np.sum(Vol_tank_list)

        Area_tank_1 = 3293         # m2
        Area_tank_2 = 0         # m2
        Area_tank_3 = 0  # m2
        Area_tank_4 = 0  # m2
        Area_tank_list = np.array([Area_tank_1, Area_tank_2, Area_tank_3, Area_tank_4])
        Area_tank = np.sum(Area_tank_list)

        Tank_info = [Type_tank, Circle_ratio, Vol_tank_list, Vol_tank, Area_tank_list, Area_tank, BOR, FR]

        Ins_thick = 300            # mm    
        Wall_thick = 12           # mm 
        M_wall =  1020000         # kg SUS weight
        Cp_wall = 0.35             # kJ/kg/C
        Cond_wall = 35             # kJ/h/m/C

        M_ins = 53000            # kg Ins weight
        Cp_ins = 1.07            # kJ/kg/C
        D_ins = 118                # kg/m3
        Cond_ins = 0.13           # kJ/h/m/C  Insulation conductivity 를 제대로 적용해야함.

        ht_NG = 420              # kJ/h/m2/C
        ht_atm = 50.0              # kJ/h/m2/C
        Tank_prop = [Ins_thick, M_wall, Cp_wall, M_ins, Cp_ins, D_ins, Cond_ins, ht_NG, ht_atm]

        sec = 0.1  # secondary barrier location in insulation

    elif Model == "SN2434": # Type B, 12K LFV

        BOR = 0.2           # %/day
        FR = 98.0            # %    
        ## Tank 정보는 Model 표 만들어서 거기서 가져오는것으로 만들기
        Vol_tank_1 = 12000         # m3
        Vol_tank_2 = 0        # m3
        Vol_tank_3 = 0    # m3
        Vol_tank_4 = 0    # m3
        Vol_tank_list = np.array([Vol_tank_1, Vol_tank_2, Vol_tank_3, Vol_tank_4])
        Vol_tank = np.sum(Vol_tank_list)

        Area_tank_1 = 3293         # m2
        Area_tank_2 = 0         # m2
        Area_tank_3 = 0  # m2
        Area_tank_4 = 0  # m2
        Area_tank_list = np.array([Area_tank_1, Area_tank_2, Area_tank_3, Area_tank_4])
        Area_tank = np.sum(Area_tank_list)

        Tank_info = [Type_tank, Circle_ratio, Vol_tank_list, Vol_tank, Area_tank_list, Area_tank, BOR, FR]

        Ins_thick = 300            # mm    
        Wall_thick = 12           # mm 
        M_wall =  1020000         # kg SUS weight
        Cp_wall = 0.35             # kJ/kg/C
        Cond_wall = 35             # kJ/h/m/C

        M_ins = 53000            # kg Ins weight
        Cp_ins = 1.07            # kJ/kg/C
        D_ins = 118                # kg/m3
        Cond_ins = 0.13           # kJ/h/m/C  Insulation conductivity 를 제대로 적용해야함.

        ht_NG = 420              # kJ/h/m2/C
        ht_atm = 50.0              # kJ/h/m2/C
        Tank_prop = [Ins_thick, M_wall, Cp_wall, M_ins, Cp_ins, D_ins, Cond_ins, ht_NG, ht_atm]

        sec = 0.1  # secondary barrier location in insulation

    elif Model == "174k_Mem":
        BOR = 0.1           # %/day
        FR = 98.5             # %

        ## Tank 정보는 Model 표 만들어서 거기서 가져오는것으로 만들기
        Vol_tank_1 = 25997         # m3
        Vol_tank_2 = 49716         # m3
        Vol_tank_3 = Vol_tank_2    # m3
        Vol_tank_4 = Vol_tank_2    # m3

        # Vol_tank_1 = 6000       # m3
        # Vol_tank_2 = 0         # m3
        # Vol_tank_3 = Vol_tank_2    # m3
        # Vol_tank_4 = Vol_tank_2    # m3

        Vol_tank_list = np.array([Vol_tank_1, Vol_tank_2, Vol_tank_3, Vol_tank_4])
        Vol_tank = np.sum(Vol_tank_list)

        Area_tank_1 = 5060         # m2
        Area_tank_2 = 7877         # m2
        Area_tank_3 = Area_tank_2  # m2
        Area_tank_4 = Area_tank_2  # m2

        # Area_tank_1 = 2300        # m2
        # Area_tank_2 = 0         # m2
        # Area_tank_3 = Area_tank_2  # m2
        # Area_tank_4 = Area_tank_2  # m2    

        Area_tank_list = np.array([Area_tank_1, Area_tank_2, Area_tank_3, Area_tank_4])
        Area_tank = np.sum(Area_tank_list)

        Tank_info = [Type_tank, Circle_ratio, Vol_tank_list, Vol_tank, Area_tank_list, Area_tank, BOR, FR]

        # return Tank_info

        Ins_thick = 270            # mm    
        Wall_thick = 1.2            # mm 

        M_wall = 383849            # kg SUS weight
        Cp_wall = 0.46             # kJ/kg/C
        Cond_wall = 35             # kJ/h/m/C

        M_ins = 1891303            # kg Ins weight
        Cp_ins = 1.071             # kJ/kg/C
        D_ins = 190                # kg/m3
        Cond_ins = 0.096           # kJ/h/m/C

        ht_NG = 37.46              # kJ/h/m2/C
        ht_atm = 16.8              # kJ/h/m2/C

        # Ins_thick = 270            # mm
        # M_wall = 31000            # kg SUS weight
        # Cp_wall = 0.46             # kJ/kg/C

        # M_ins = 170000            # kg Ins weight
        # Cp_ins = 1.07             # kJ/kg/C
        # D_ins = 118               # kg/m3
        # Cond_ins = 0.08           # kJ/h/m/C

        # ht_NG = 37.46              # kJ/h/m2/C
        # ht_atm = 16.8              # kJ/h/m2/C

        Tank_prop = [Ins_thick, M_wall, Cp_wall, M_ins, Cp_ins, D_ins, Cond_ins, ht_NG, ht_atm]

        sec = 0.1  # secondary barrier location in insulation
        # return  Tank_prop

        # def TypeC_7_5k() : # SN2430

        #     BOR = 0.2           # %/day
        #     FR = 90             # %    
        #     ## Tank 정보는 Model 표 만들어서 거기서 가져오는것으로 만들기
        #     Vol_tank_1 = 3750         # m3
        #     Vol_tank_2 = 3750         # m3
        #     Vol_tank_3 = 0    # m3
        #     Vol_tank_4 = 0    # m3
        #     Vol_tank_list = np.array([Vol_tank_1, Vol_tank_2, Vol_tank_3, Vol_tank_4])
        #     Vol_tank = np.sum(Vol_tank_list)

        #     Area_tank_1 = 2874/2         # m2
        #     Area_tank_2 = 2874/2         # m2
        #     Area_tank_3 = 0  # m2
        #     Area_tank_4 = 0  # m2
        #     Area_tank_list = np.array([Area_tank_1, Area_tank_2, Area_tank_3, Area_tank_4])
        #     Area_tank = np.sum(Area_tank_list)

        #     Tank_info = [Type_tank, Circle_ratio, Vol_tank_list, Vol_tank, Area_tank_list, Area_tank, BOR, FR]

        #     Ins_thick = 300            # mm    
        #     M_wall = 550000            # kg SUS weight
        #     Cp_wall = 0.35             # kJ/kg/C

        #     M_ins = 44000            # kg Ins weight
        #     Cp_ins = 1.071             # kJ/kg/C
        #     D_ins = 118                # kg/m3
        #     Cond_ins = 0.08           # kJ/h/m/C

        #     ht_NG = 37.46              # kJ/h/m2/C
        #     ht_atm = 19.5              # kJ/h/m2/C
        #     Tank_prop = [Ins_thick, M_wall, Cp_wall, M_ins, Cp_ins, D_ins, Cond_ins, ht_NG, ht_atm]
    U = 1/(1/ht_NG + Wall_thick/1000/Cond_wall + Ins_thick/1000/Cond_ins + 1/ht_atm)
    print(Model, Vol_tank, Ins_thick)

[JungHulk]

# -*- coding: utf-8 -*-

"""
22. 06. 27
"""

import os, numpy as np
import math
import imageio
import matplotlib.pyplot as plt
from win32com.client import Dispatch
import pandas as pd
import numpy as np
# from scipy import integrate
import pickle
import csv
import openpyxl
import sys
from ctREFPROP.ctREFPROP import REFPROPFunctionLibrary
from scipy import integrate
import Tank_information as T
import time
Tank = T.Cargo_tank()
seconds = 3
while seconds > 0:
    print("Check Model :", seconds, "sec")
    time.sleep(0.5)
    seconds = seconds - 1

os.environ['RPPREFIX'] = r'E:\REFPROP'

RP = REFPROPFunctionLibrary(os.environ['RPPREFIX'])
RP.SETPATHdll(os.environ['RPPREFIX'])

SI = RP.GETENUMdll(0,"SI WITH C").iEnum
# Property = "P;T;D;H;M;CP;Phase;Qmass"
# Prop_list = "P;T;D;H;M;CP;Phase;Qmass;VIS;TCX;PRANDTL;KV;Heatvapz_P"    
# Property_total = "T;P;D;H;CP;Phase;Qmass;"
Prop_list = "P;T;D;H;M;CP;Qmass;VIS;TCX;PRANDTL;KV;Heatvapz_P;Phase;"     # VISCOSITY 는 100만 나누어야 함.
Prop_index = np.array(Prop_list.replace(';',' ').split())

# df1 = pd.read_pickle('Pipetable.p')    

def Heat_transfer(Supply_condition, Init_cond, P_mode, T_ins, T_ins_ave, Ins_node, dt):

    #################### Initial Calculation condition
    # Init_cond = Init_condition(Comp_name, Comp_init, P_init, Surface, Filling, Vol_liq, Vol_vap, Step, LastStep, Target)
    T_vap = Init_cond[2][1]
    h = Init_cond[2][3]
    d = Init_cond[2][2]
    Phase = Init_cond[5]
    Qmass = Init_cond[2][6]

    h_liq = Init_cond[3][3]
    d_liq = Init_cond[3][2]

    h_vap = Init_cond[4][3]
    d_vap = Init_cond[4][2]

    Comp = Comp_init
    Comp_liq = Init_cond[0]
    Comp_vap = Init_cond[1]

    # P = RP.REFPROPdll(Comp_name, "DH" , "P", SI, 1, 0, d, h, Comp).Output[0]
    # print('Pressure : ', P)
    Supply = Supply_condition
    val = Supply[4].find("Warm") #Operation check Warm-up 시에는 열공급 효율 70%로 진행
    val2 = Supply[4].find("Inert") #Operation check Inert 시에는 O2, CH4 composition print
    val_cool = Supply[4].find("Cool") #Operation check CD 시에는 spray 누적량 표기, Exhaust 량 표기

    dt = dt
    # dtt = 2 * dt

    T_wall = T_vap
    Mass_tank = Volume * Init_cond[2][2]
    Mass_vap = Mass_tank * Qmass
    Mass_liq = Mass_tank - Mass_vap

    T_ins_var = 0
    T_ins = T_ins
    Ins_node = Ins_node
    T_ins_temp = np.zeros(len(T_ins))

    Mass_exhaust_init = 0 # 100 kg/h
    Vol_exhaust_init = 0

    Amount_spray = 0 # kg/h accumulate

    T_ins_sec = (T_ins[loc_sec] + T_ins[loc_sec_1]) / 2

    Q_tank = U * Surface * (T_atm - T_vap)

    ################### Data save #############################################
    strFormat = '%'

    T_ins_table = np.append(T_ins, T_ins_ave) 
    T_ins_table = np.append(T_ins_table, T_ins_sec)
    Prop_table = np.array([0, P_init, T_vap, T_wall, d, h, Mass_tank, Phase, Qmass, Mass_exhaust_init, Vol_exhaust_init, Amount_spray])

    print(Comp_init)
    print(Comp_liq)
    print(Comp_vap)
    Comp_table = np.array([Comp_init, Comp_liq, Comp_vap])
    heatleak = []

    for i in range (Calculation):

        # if i < 100:  # 초기 계산에만 time step 줄여서 계산하다가 100번 iteration 후에는 time step 10배로
        #     dt = dt
        # else:
        #     dt = dtt
        # if Operation == "Cool-down":
        #     if T_vap < -140:
        #         break
        # elif Operation == "Cool-down":

        if val == -1:
            dh = (-1 * (Supply[2][i] * (h - Supply[1][3]) + ht_NG * Surface * (T_vap - T_wall)) * dt) / Mass_tank

        else:  #Operation check :  Warm-up 시에는 열공급 효율 70%로 진행
            dh = (-1 * (Supply[2][i]*0.7 * (h - Supply[1][3]) + ht_NG * Surface * (T_vap - T_wall)) * dt) / Mass_tank
        # print(i*dt+dt, T_vap, dh, Mass_tank)

        dT_wall = ((ht_NG * Surface * (T_vap - T_wall) + ht_atm * Surface * (T_atm - T_ins[-1])) * dt - M_ins * Cp_ins* T_ins_var) / (M_wall * Cp_wall)
        Comp = list(map(lambda x, y: (x * Mass_tank + y * Supply[2][i] * dt) / (Mass_tank + Supply[2][i]*dt), Comp, Supply[0])) 

        h_pre = h # 이전 스텝 h값
        # 압력항에 대한 것도 같이 넣어줘야 할것 같음. 
        if P_mode == 0: # 0: Const, 1: Various
            h = h + dh
            d = RP.REFPROPdll(Comp_name, "PH" , "D", SI, 1, 0, P_init, h, Comp).Output[0]
            Mass_exhaust = (Mass_tank + (Supply[2][i] * dt) - (Volume * d)) / dt
            Vol_exhaust = Mass_exhaust / d
            ### Vol_exhause = (Mass_tank/d - Mass_tank/d_pre + Spray_mass * dt/d) / dt
            Mass_tank = Volume * d
            P = P_init

            T_vap = np.array(RP.REFPROPdll(Comp_name,"DH","T",SI,1,1,d,h,Comp).Output[0])
            Qmass = np.array(RP.REFPROPdll(Comp_name,"DH","Qmass",SI,1,1,d,h,Comp).Output[0])
            Phase = np.array(RP.REFPROPdll(Comp_name,"DH","Phase",SI,1,1,d,h,Comp).hUnits)
            T_wall = T_wall + dT_wall  

        else:
            Mass_exhaust = 0   ## user가 수정

            h = ((-1 * (Supply[2][i] * (- Supply[1][3]) + ht_NG * Surface * (T_vap - T_wall)) * dt) + Mass_liq * h_liq + (Mass_vap - Mass_exhaust * dt) * h_vap) / Mass_tank
            # print(((ht_atm * Surface) * (T_atm - T_ins[-1]))/509)

            Mass_tank = Mass_tank + Supply[2][i]*dt-Mass_exhaust
            d = Mass_tank / Volume
            P = RP.REFPROPdll(Comp_name, "DH" , "P", SI, 1, 1, d, h, Comp).Output[0]
            Vol_exhaust  = Mass_exhaust / d   ## user가 수정

            T_vap = RP.REFPROPdll(Comp_name,"DH","T",SI,1,1,d,h,Comp).Output[0]
            Qmass = RP.REFPROPdll(Comp_name,"DH","Qmass",SI,1,1,d,h,Comp).Output[0]
            Phase = np.array(RP.REFPROPdll(Comp_name,"DH","Phase",SI,1,1,d,h,Comp).hUnits)
            T_wall = T_wall + dT_wall       

            Comp_liq = np.array(RP.REFPROPdll(Comp_name,"DT","xmassliq",SI,1,0,d,T_vap,Comp).Output[0:Comp_name.count(';')])
            Comp_vap = np.array(RP.REFPROPdll(Comp_name,"DT","xmassvap",SI,1,0,d,T_vap,Comp).Output[0:Comp_name.count(';')])

            Mass_vap = Mass_tank * Qmass
            Mass_liq = Mass_tank - Mass_vap

            if Qmass >0 and Qmass < 1: # two-phase

                P = RP.REFPROPdll(Comp_name, "DH" , "P", SI, 1, 1, d, h, Comp).Output[0]

                # d_vap    = RP.REFPROPdll(Comp_name,"PQmass","Dvap",SI,1,1,P, Qmass,Comp_vap).Output[0]
                # d_liq    = RP.REFPROPdll(Comp_name,"PQmass","Dliq",SI,1,1,P, Qmass,Comp_liq).Output[0]
                # h_vap    = RP.REFPROPdll(Comp_name,"PD","H",SI,1,1,P,d_vap, Comp_vap).Output[0]
                # h_liq    = RP.REFPROPdll(Comp_name,"PD","H",SI,1,1,P,d_liq, Comp_liq).Output[0]

                d_vap    = RP.REFPROPdll(Comp_name,"Ph","Dvap",SI,1,1,P, h, Comp).Output[0]
                d_liq    = RP.REFPROPdll(Comp_name,"Ph","Dliq",SI,1,1,P, h, Comp).Output[0]
                # h_vap    = RP.REFPROPdll(Comp_name,"PD","H",SI,1,1,P,d_vap, Comp_vap).Output[0]
                # h_liq    = RP.REFPROPdll(Comp_name,"PD","H",SI,1,1,P,d_liq, Comp_liq).Output[0]
                h_vap    = RP.REFPROPdll(Comp_name,"PD","Hvap",SI,1,1,P,d, Comp).Output[0]
                h_liq    = RP.REFPROPdll(Comp_name,"PD","Hliq",SI,1,1,P,d, Comp).Output[0]               

            else:
                d_vap    = np.array(RP.REFPROPdll(Comp_name,"PT","Dvap",SI,1,1,P,T_vap,Comp_vap).Output[0])
                d_liq    = np.array(RP.REFPROPdll(Comp_name,"PT","Dliq",SI,1,1,P,T_vap,Comp_liq).Output[0])
                h_vap    = np.array(RP.REFPROPdll(Comp_name,"PD","h",SI,1,1,P,d_vap,Comp_vap).Output[0])
                h_liq    = np.array(RP.REFPROPdll(Comp_name,"PD","h",SI,1,1,P,d_liq,Comp_liq).Output[0])

        # Prop = Properties(Comp_name, Comp, "DH", d, h)[0]
        # Prop = np.array(RP.REFPROPdll(Comp_name,"DH",Prop_list,SI,1,0,d,h,Comp).Output[0:Prop_list.count(';')])

        Comp_array = np.array([Comp, Comp_liq, Comp_vap])
        Comp_table = np.vstack([Comp_table, Comp_array])        
        # Comp_table = np.vstack([Comp_table, Comp]) 
        for j in range (len(T_ins)-2):
            k = j+1
            T_ins_temp[0] = T_wall
            T_ins_temp[k] = diffusivity * dt / Ins_interval**2 * (T_ins[j] + T_ins[j+2]) + (1- 2* (diffusivity * dt / Ins_interval**2)) * T_ins[j+1]
            T_ins_temp[len(T_ins)-1] = Cond_ins/ ht_atm * (T_ins[len(T_ins)-2] - T_ins[len(T_ins)-1]) / (Ins_node[len(T_ins)-1] - Ins_node[len(T_ins)-2]) + T_atm

        T_ins_sec = (T_ins[loc_sec] + T_ins[loc_sec_1]) / 2
        T_ins_temp_ave = simpson(T_ins_temp,len(T_ins)-1)
        T_ins_temp_table = np.append(T_ins_temp, T_ins_temp_ave)
        T_ins_temp_table = np.append(T_ins_temp_table, T_ins_sec)
        T_ins_table = np.vstack([T_ins_table, T_ins_temp_table]) # Time step 별 Insul node 별 테이블화

        Amount_spray = Amount_spray + Supply[2][i]*dt
        Prop_tank = np.array([np.round(i*dt+dt,2), np.round(P,4), np.round(T_vap,4), np.round(T_wall,4), np.round(d,2), np.round(h,4), np.round(Mass_tank,2), Phase, np.round(Qmass,4), np.round(Mass_exhaust,2), np.round(Vol_exhaust,2), Amount_spray])
        # print (Prop_tank)
        Prop_table = np.vstack([Prop_table, Prop_tank])

        T_ins_var = T_ins_temp_ave - T_ins_ave
        # print(T_ins_temp_ave, T_ins_ave, T_ins_var)
        T_ins = T_ins_temp
        T_ins_ave = T_ins_temp_ave
        print('Time :', "%-8s" % np.round(i*dt+dt,2), 'Press. :', "%-8s" % np.round(P,6), 'T_vap :', "%-8s" % np.round(T_vap,6), 'T_wall :', "%-8s" % np.round(T_wall,2), 'T_sec :', np.round(T_ins_sec,2))

        if val_cool is not -1:
            print('Amount of Spray  :  ', "%-10s" % np.round(Amount_spray,2), "Exhaust  :  ", "%-10s" % np.round(Mass_exhaust,2))
        ### dew point 
        ### https://en.wikipedia.org/wiki/Dew_point
        ### https://en.wikipedia.org/wiki/Arden_Buck_equation
        if val2 == -1:
            dew_point = 0
        else:
            water_mole = RP.REFPROPdll(Comp_name,"DT","Xmole",SI,1,1,d,T_vap,Comp).Output[10]
            print(water_mole)
            RH = (P*water_mole)/RP.REFPROPdll("Water","TQmass","P",SI,1,1,T_vap,0,[1.0]).Output[0]  ## 상대 습도
            print(RH)
            dew_point = 243.04*(np.log(RH)+(17.625*T_vap/(243.04+T_vap)))/(17.625-np.log(RH)-(17.625*T_vap/(243.04+T_vap)))

            print("Methane : " , "%-10s"%np.round(Comp[1]*100,4), "Oxygen : ", "%-10s"%np.round(Comp[8]*100,4), "N2 :" , Amount_spray/1000)    
            print("Dew Point : ", dew_point)
        Heat_to_gas = Mass_tank * (h - h_pre) / dt # kJ/h
        Heat_to_pri = M_wall * Cp_wall * T_ins_var / dt   # kJ/h

        ################################ Plot ###########################################
        ### 시간에 따라 지속적으로 plot
        # fig = plt.figure(figsize=(10,10))
        # ax = fig.add_subplot()

        # x_axis.append(i * dt)
        # y_axis.append(T_vap)
        # y_axis2.append(T_wall)
        # y_axis3.append(T_ins_ave)
        # # ax.plot(x_axis,y_axis, color = 'b', linewidth = 4.0)
        # # ax.plot(x_axis,y_axis2, color = 'r', linewidth = 4.0)
        # # ax.plot(x_axis,y_axis3, color = 'green', linewidth = 4.0)
        # plt.plot(x_axis, y_axis, 'b', x_axis, y_axis2, 'r', x_axis, y_axis3, 'green', linewidth = 4.0 )

        # ax.patch.set_facecolor('#c6d9f1') # axes 배경색
        # # plt.setp(line, color = 'b', linewidth = 4.0)
        # plt.xlabel('Time (hr)', size=15)
        # plt.ylabel('Temperature (C)', size=15)
        # # plt.rc('axes', labelsize = 70)

        # plt.xlim(0,)
        # plt.ylim(-180,40)
        # plt.grid(True)
        # plotting = int(i % 100)        

        ### 100 step 마다 plot
        # plt.plot(1,1)
        # plt.figure(figsize=(10,10))
        plt.rcParams["figure.figsize"] = (10,10)
        ax = plt.gca()
        ax.set_facecolor('#c6d9f1') # axes 배경색
        x_axis.append(i * dt)
        y_axis.append(T_vap)
        y_axis2.append(T_wall)
        y_axis3.append(T_ins_ave)
        y_axis4.append(T_ins_sec)
        # plt.subplot(211)
        plt.plot(x_axis, y_axis, 'b', x_axis, y_axis2, 'r', x_axis, y_axis3, 'green', x_axis, y_axis4, 'Orange', linewidth = 4.0 )

        plt.xlabel('Time (hr)', size=15)
        plt.ylabel('Temperature (C)', size=15)

        plt.xlim(0, Calculation*dt)
        plt.ylim(-180,50)
        plt.grid(True)
        plotting = int(i % 100)
        gif = int(i % 10)
        ### GIF 만들기
        # if gif == 0:
        #     plotname = f'Plot{i}.png'
        #     plotnames.append(plotname)
        #     plt.savefig(plotname, bbox_inches = 'tight', pad_inches = 0)

        # if i == (Calculation-1):
        #     name = f'Last.jpeg'
        #     plt.savefig(name, bbox_inches = 'tight', pad_inches = 0 )

        if  plotting == 0:
            plt.show()

        ### Insulation 온도 분포 (누적, 또는 나중에 Table 해서 그냥 한번만 뽑기)
        # if  plotting == 0:
        #     x2_axis = Ins_node
        #     y2_axis = T_ins
        #     plt.subplot(212)

        #     plt.rcParams["figure.figsize"] = (10,10)
        #     ax = plt.gca()
        #     ax.set_facecolor('#c6d9f1') # axes 배경색
        #     plt.xlim(0,Ins_thick/1000)
        #     plt.ylim(-170,50)
        #     plt.xlabel('Insulation Thickness (m)', size=15)
        #     plt.ylabel('Temperature (C)', size=15)
        #     plt.grid(True)

        #     print(x2_axis, y2_axis)
        #     plt.plot(x2_axis, y2_axis, linewidth = 4.0)

        #     plt.show()

    return Prop_table, T_ins_table, Comp_table

    # Exhaust vapor 계산하기
    # Mass_conv = Mass_supply * dt + Mass_tank - Mass_exhaust

    # for i in (len(T_ins)-1):
    #     k = i+1
    #     T_ins_temp[i] = T_wall
    #     T_ins_temp[k] = diffusivity * dt / Ins_interval**2 * (T_ins[i] + T_ins[i+2]) + (1- 2* (diffusivity * dt / Ins_interval**2)) * T_ins[i+1]  

# for i in range (0,12/3600+t_step,t_step):

def Supply_condition(Comp_name, Mode_operation, Input_prop, Input1, Input2, Mass_supply, Calculation):
# 0 'Inerting', 1'Inerting_N', 2'Gassing-up', 3'Cool-down', 4'Cool-down_N', 5'Warm-up', 6'Warm-up_N', 7'Aeration', 8'Hold']    
    if Mode_operation == Mode_list[0]:
        Comp_supply = np.array([0.85, 0, 0, 0, 0, 0, 0, 0, 0.01, 0.14, 0])  

    elif Mode_operation == Mode_list[1] or Mode_operation == Mode_list[4] or Mode_operation == Mode_list[6]:
        Comp_supply = np.array([1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

    elif Mode_operation == Mode_list[2] or Mode_operation == Mode_list[3] or Mode_operation == Mode_list[5] or Mode_operation == Mode_list[8]:
        Comp_supply = Comp_init

    else: # Aeration
        Comp_supply = np.array([0.78, 0, 0, 0, 0, 0, 0, 0, 0.21, 0.01, 0]) 

    Value = Properties(Comp_name, Comp_supply, Input_prop, Input1, Input2)[0]
    Volume_supply = Mass_supply / Value[2]

    ## Supply ramp-up
    Cal_time = Calculation
    Supply_mass = np.full(Cal_time, Mass_supply)

    for i in range (11):
        if i < 6:
            Supply_mass[i] = 1/(10-i) * Mass_supply
        else:
            Supply_mass[i] = (i-5) / 5 * Mass_supply

    return Comp_supply, Value, Supply_mass, Volume_supply, Mode_operation

def Properties(Comp_name, Comp_value, Input_prop, input1, input2):
    Prop = []
    Value = []

    for k in Prop_index:   # P;T;D;H;M;CP;Qmass;VIS;TCX;PRANDTL;KV;Heatvapz_P;Phase;
        r = RP.REFPROPdll(Comp_name, Input_prop, k, SI, 1, 0, input1, input2, Comp_value)
        Prop.append([k,r.hUnits,r.Output[0]])
        Value.append(r.Output[0])

    DY = Value[2] * Value[11] / (100**2) # Dynamic viscousity

    df_prop = pd.DataFrame(Prop, index = Prop_index, columns=['Prop','Units','Values'])
    df_prop = df_prop.drop(['Prop'], axis=1)
    df_prop.loc['DV'] = ['kg/m-s', DY]

    return Value, df_prop # Value, dataframe 

def Init_condition(Comp_name, Comp_init, P_init, Filling, Step, LastStep, Target):

    Qmass = 0.01  # Initial Qmass
    direction = 0 # Try and Error Initial direction

    D_tank = RP.REFPROPdll(Comp_name,"PQmass","D",SI,1,0,P_init,Qmass,Comp_init).Output[0]
    D_liq = RP.REFPROPdll(Comp_name,"PQmass","Dliq",SI,1,0,P_init,Qmass,Comp_init).Output[0]
    D_vap = RP.REFPROPdll(Comp_name,"PQmass","Dvap",SI,1,0,P_init,Qmass,Comp_init).Output[0]
    Vol_liq = Volume * Filling/100
    Vol_vap = Volume - Vol_liq

    Mass_tank = Volume * D_tank
    Mass_liq = Vol_liq * D_liq
    Mass_vap = Vol_vap * D_vap

    Mass_conv = Mass_tank - Mass_liq - Mass_vap # 질량보존 확인항

    if Mass_conv < Target:
        direction = 1
    else:
        direction = -1
    if  Mass_conv == Target:
        sys.exit(0) # sys.exit(0) 프로그램 정상 종료, (1) 강제 종료

    if  Filling == 0:

        Comp_liq_mass = np.array(RP.REFPROPdll(Comp_name,"PT","xmassliq",SI,1,0,P_init,T_set_vapor,Comp_init).Output[0:Comp_name.count(';')])
        Comp_vap_mass = np.array(RP.REFPROPdll(Comp_name,"PT","xmassvap",SI,1,0,P_init,T_set_vapor,Comp_init).Output[0:Comp_name.count(';')])

        Pro_init_mass = np.array(Properties(Comp_name, Comp_init, "PT", P_init, T_set_vapor)[0])
        Pro_liq_mass = Pro_init_mass
        Pro_vap_mass = Pro_init_mass       

        Phase = RP.REFPROPdll(Comp_name,"PT","Phase",SI,1,0,P_init,T_set_vapor,Comp_init).hUnits

    else:
        for i in range(1,10000,1):

            Qmass = Qmass + Step * direction  

            if i % 100 == 0: ## 수렴 가속화
                print('=====================================================')
                print('Qmass reset')
                print('=====================================================')
                Qmass = Mass_vap / Mass_tank
                # Step = Step * 10

            if Mass_conv < Target:
            ### Qmass 가 음수가 되는 경우 피하기 위한 함수
                if Qmass < 0:
                   if direction == -1:
                       if Step > LastStep:
                            # Step = Step / 10
                            direction = direction * -1
            ##############################################

                elif direction == 1:
                    if Step > LastStep:
                        Step = Step / 10

                        direction = direction * -1          

            else:
                if direction == -1:
                    if Step > LastStep:
                        Step = Step / 10
                        direction = direction * -1

            if Step < LastStep:

                break
            aa = Mass_vap / Mass_tank
            # print(Step, aa, direction)
            # print(D_tank, D_liq, D_vap)
            # print(i, Qmass, Step)
            print("Calculating...  i: {} Mass_conv {} Tank {} Liq {} Vap {} Qmass: {}".format(i,round(Mass_conv,8), round(Mass_tank,2), round(Mass_liq,2), round(Mass_vap,2), round(Qmass,15)))
            D_tank = RP.REFPROPdll(Comp_name,"PQmass","D",SI,1,0,P_init,Qmass,Comp_init).Output[0]
            D_liq = RP.REFPROPdll(Comp_name,"PQmass","Dliq",SI,1,0,P_init,Qmass,Comp_init).Output[0]
            D_vap = RP.REFPROPdll(Comp_name,"PQmass","Dvap",SI,1,0,P_init,Qmass,Comp_init).Output[0]

            Mass_tank = Volume * D_tank
            Mass_liq = Vol_liq * D_liq
            Mass_vap = Vol_vap * D_vap

            Mass_conv = Mass_tank-Mass_liq-Mass_vap
        print('Done', i)
        print('Qmass  : ', Qmass)
        print('Mass_conv  : ', Mass_conv)

        Comp_liq_mass = np.array(RP.REFPROPdll(Comp_name,"PQmass","xmassliq",SI,1,0,P_init,Qmass,Comp_init).Output[0:Comp_name.count(';')])
        Comp_vap_mass = np.array(RP.REFPROPdll(Comp_name,"PQmass","xmassvap",SI,1,0,P_init,Qmass,Comp_init).Output[0:Comp_name.count(';')])

        Pro_init_mass = np.array(Properties(Comp_name, Comp_init, "PQmass", P_init, Qmass)[0])
        Pro_liq_mass = np.array(Properties(Comp_name, Comp_liq_mass, "PQmass", P_init, 0)[0])
        Pro_vap_mass = np.array(Properties(Comp_name, Comp_vap_mass, "PQmass", P_init, 1)[0])

        Phase = RP.REFPROPdll(Comp_name,"PQmass","Phase",SI,1,0,P_init,Qmass,Comp_init).hUnits

    return Comp_liq_mass, Comp_vap_mass, Pro_init_mass, Pro_liq_mass, Pro_vap_mass, Phase

def Unitconv_nominal(Comp_name, Comp, P, T, iflag, Value): # iflag 0: Nominal flowrate, 1: Volumetric flowrate, 2 : Mass flowrate

    d = RP.REFPROPdll(Comp_name, "PT" , "D", SI, 1, 0, P, T, Comp).Output[0]
    d_no = RP.REFPROPdll(Comp_name, "PT" , "D", SI, 1, 0, 0.101325, 0, Comp).Output[0]

    if iflag == 0:
        Nominal = Value
        Mass    = Nominal * d_no
        Vol     = Mass / d

    elif iflag == 1:
        Vol     = Value
        Mass    = Vol * d
        Nominal = Mass / d_no

    elif iflag == 2:
        Mass    = Value
        Vol     = Mass / d
        Nominal = Mass / d_no

    print("Nominal : ", np.round(Nominal,2), "Nm3/h")
    print("Volume flow : ", np.round(Vol,2), "m3/h")
    print("Mass flow : ", np.round(Mass,2), "kg/h")
    return Nominal, Vol, Mass

def BOG_calculation(Tank_volume, Filling_Ratio, BOR, Comp_name, Comp_LNG): # methane, LNG, Revised_BOR
    Vol = Tank_volume
    d_methane = RP.REFPROPdll("Methane", "PQmass", "D", SI, 1, 0, P_init, 0, [1.0]).Output[0]
    LH_methane = RP.REFPROPdll("Methane", "PQmass", "Heatvapz_P", SI, 1, 0, P_init, 0, [1.0]).Output[0]

    d = RP.REFPROPdll(Comp_name, "PQmass", "D", SI, 1, 0, P_init, 0, Comp_LNG).Output[0]

    ## LNG 기화된 Vap 의 LH 사용
    Comp_vap_LNG = RP.REFPROPdll(Comp_name,"PQmass","XmassVAP",SI,1,0,P_init,0.01,Comp_LNG).Output[0:len(Comp_LNG)]
    LH = RP.REFPROPdll(Comp_name, "PQmass", "Heatvapz_P", SI, 1, 0, P_init, 0, Comp_vap_LNG).Output[0]

    # print(Comp_vap_LNG)
    # print(LH)
    BOG_methane = Vol * FR/100 * d * BOR /24 /100
    BOG_LNG     = BOG_methane * LH_methane * d_methane / LH / d

    Revised_BOR = BOR * LH_methane * d_methane / LH / d

    print("BOG Methane :", "%-5s" % round(BOG_methane,2), "%-5s" % "kg/h")
    print("BOG_LNG :", "%-5s" % round(BOG_LNG,2), "kg/h")
    return BOG_methane, BOG_LNG, Revised_BOR       # 메탄 BOG(kg/h) LNG BOG(kg/h), IMO 수정된 BOR

# def BOR_check(BOG, Comp_name, Comp_init, P_init, Filling, Step, LastStep, Target, Mode_list, Sim_time):

#     Init_cond = Init_condition(Comp_name, Comp_init, P_init, Filling, Step, LastStep, Target)

#     Operation = Mode_list[6]
#     Supply = Supply_condition(Comp_name, Operation, "PT", 0.11, 60, 0, Sim_time)

#     Hold = Heat_transfer(Supply, Init_cond, 0, T_ins, T_ins_ave, Ins_node, dt)
#     Exhaust = Hold[0][9]

    # if Exhaust < BOG:
    #     direction = 1
    # else:
    #     direction = -1
    # if  round(Exhaust,1) < round(BOG,1)*1.1 and round(Exhaust,1) > round(BOG,1)*0.9:
    #     sys.exit(0) # sys.exit(0) 프로그램 정상 종료, (1) 강제 종료

    # for i in range(1,10000,1):
    #     print(ht_atm, Exhaust, BOG)
    #     ht_atm = ht_atm + Step * direction  
    #     Exhaust = Heat_transfer(Supply, Init_cond, 0, T_ins, T_ins_ave, Ins_node, dt)[0][9]
    #     if Exhaust < BOG:
    #     ### Qmass 가 음수가 되는 경우 피하기 위한 함수
    #         if ht_atm < 0:
    #            if direction == -1:
    #                if Step > LastStep:
    #                    Step = Step / 10
    #                    direction = 1
    #     ##############################################

    #         elif direction == 1:
    #             if Step > LastStep:
    #                 Step = Step / 10

    #                 direction = direction * -1          

    #     else:
    #         if direction == -1:
    #             if Step > LastStep:
    #                 Step = Step / 10
    #                 direction = direction * -1

    #     if Step < LastStep:
    #         break

    # return Exhaust

'''''''''''''''''''''''''''
        Initial Condition
'''''''''''''''''''''''''''
Comp_name =          "Nitrogen;Methane;Ethane;Propane;Butane;Isobutane;Pentane;Isopentane;Oxygen;CO2;Water;"  
Comp_name_index =   np.array(Comp_name.replace(';',' ').split())
Comp_init = np.array([1.0,         0.0,    0,    0,       0,        0,     0,        0,       0,    0,   0]) # pure N2
# Comp_init = np.array([0.0,         1,    0,    0,       0,        0,     0,        0,       0,    0,   0])
# Pure_methane  = np.array([0.0,         1,    0,    0,       0,        0,     0,        0,       0,    0,   0])
# Comp_init = np.array([0.0079, 0.9205, 0.0547, 0.0133, 0.0037, 0, 0, 0, 0, 0]) # Mass base
# Comp_init = np.array([0.751,         0.0,    0,    0,       0,        0,     0,        0,       0.2319,    0.0152,   0.0019]) # Air

T_atm = 25 # C  
P_init = 0.126325 # Mpa
T_set_vapor = -160

Volume =  Tank.Vol_tank # m3
Surface = Tank.Area_tank # m2
Filling = 0 # %

FR = Tank.FR
BOR = Tank.BOR
'''''''''''''''''''''''''''''
      Insulation information 
'''''''''''''''''''''''''''''
M_wall = Tank.M_wall # kg
Cp_wall = Tank.Cp_wall # kJ/kg/C
Cond_wall = Tank.Cond_wall # kJ/h/m/C
Wall_thick = Tank.Wall_thick # mm

M_ins = Tank.M_ins # kg
Cp_ins = Tank.Cp_ins # kJ/kg/C
D_ins = Tank.D_ins # kg/m3

Cond_ins = Tank.Cond_ins # kJ/h/m/C
Ins_thick = Tank.Ins_thick # mm

ht_NG = Tank.ht_NG      # kJ/h/m2/C
ht_atm = Tank.ht_atm    # kJ/h/m2/C

U = Tank.U

'''''''''''''''''''''''''''
         Equilibrium Mode
'''''''''''''''''''''''''''
# Qmass = 0.1 # initial value funcion 에 기입
Step = 0.2 
LastStep = 0.0000000000001
Target = 0.001
# direction = 0  

'''''''''''''''''''''''''''
         Plot
'''''''''''''''''''''''''''
x_axis = []
y_axis = []
y_axis2 = []
y_axis3 = []
y_axis4 = [] # sec barrier

x2_axis = []
y2_axis = []

filenames = []
plotnames = []

##############################################################
####################### Definition ###########################
##############################################################

# Input_prop = 'PT', 'PH' 등등
Init_cond = Init_condition(Comp_name, Comp_init, P_init, Filling, Step, LastStep, Target)

T_vap = Init_cond[4][1]
T_lng = Init_cond[2][1]
T_liq = Init_cond[3][1]
T_wall = T_vap

Mass_tank = Volume * Init_cond[2][2]
Mass_vap = Mass_tank * Init_cond[2][6]
Mass_liq = Mass_tank - Mass_vap

'''''''''''''''''''''''''''
Select Mode
''''''''''''''''''''''''''' 
Mode_list = ['Inerting', 'Inerting_N', 'Gassing-up', 'Cool-down', 'Cool-down_N', 'Warm-up', 'Warm-up_N', 'Aeration', 'Hold']
Mode_pressure = "Constant" # Constant, Variable

'''''''''''''''''''''''''''''''''''''''
 Insulation Initial condition
'''''''''''''''''''''''''''''''''''''''
n = 12  # Insulation node 수
Ins_interval = Ins_thick/1000 / n # Node 별 거리
Ins_node = np.arange(0,Ins_thick/1000+Ins_interval/100,Ins_interval) # Insulation 각 위치
T2 = (T_wall + (ht_atm * Ins_thick/1000/Cond_ins)*T_atm) / (1+(ht_atm*Ins_thick /1000/Cond_ins))
#Ins Outer temp.
T_ins_step = (T2 - T_wall) / (Ins_thick/1000/Ins_interval) # ins del_T
T_ins = np.zeros(n+1)

### Insulation secondary barrier 위치 정하기
Ins_sec = Tank.sec
if np.any(Ins_node == Ins_sec) == True:    # Insulation node 와 secondary barrier 가 같은경우,
    loc_sec   = np.where(Ins_node == Ins_sec)[0][0]  #Insulation location
    loc_sec_1 = np.where(Ins_node == Ins_sec)[0][0]
    # T_ins_sec = (T_ins[loc_sec] + T_ins[loc_sec_1]) / 2
else:                                      # 다를 경우는, 근처 2개 값의 평균값을 사용
    nearest = Ins_node.flat[np.abs(Ins_node - Ins_sec).argmin()]
    if nearest > Ins_sec:
        Ins_sec_1 = Ins_sec - 0.01
    elif nearest < Tank.sec:
        Ins_sec_1 = Ins_sec + 0.01
    nearest_1 = Ins_node.flat[np.abs(Ins_node - Ins_sec_1).argmin()]
    loc_sec = np.where(Ins_node == nearest)[0][0]
    loc_sec_1 = np.where(Ins_node == nearest_1)[0][0]
    # T_ins_sec = (T_ins[loc_sec] + T_ins[loc_sec_1]) / 2

# Insulation 초기온도 조건
if T_ins_step == 0:
    T_ins = np.full(n+1, T_wall)
else:
    T_ins = np.linspace(T_wall, T2, n+1)

def simpson(array, n): # simpson's rule
    integration = array[0] + array[-1]

    for i in range(1,n):
        if i%2 ==0:
            integration = integration + 2*array[i]
        else:
            integration = integration + 4*array[i]
    integration = integration / (len(array)-1) /3

    return integration

T_ins_ave = simpson(T_ins,n)

dt_ins = 0

diffusivity = Cond_ins / (Cp_ins * D_ins)

## 제어기? 혹은 기준설정 후 자동으로 종료되게 만들기
# def Operation_criteria(Mode_list, Operation, Temperature, Composition):
# ### [0] Inerting 
# ### [1] Inerting_N
# ### [2] Gassing-up
# ### [3] Cool-down
# ### [4] Cool-down_N
# ### [5] Warm-up
# ### [6] Warm-up_N
# ### [7] Aeration
#     Operation = Mode_list
#     if Operation.find("Cool") is not -1:
#         if Temperature < -140:
#             break

#     elif Operation.find("Inert") is not -1:
#         if Composition < 0.02:
#             break

# Prop = Properties(Comp_name, Comp_init, "PT", 0.12, 0)
"""""""""""""""""""""""
Calculation
"""""""""""""""""""""""
# Supply 조건 지정.
# Input_prop = "PQmass" # 2가지 고르기 (input에서)
# Input1 = 0.106 # MPaA  Input 1 값
# Input2 = 0   # Input 2 값
simulation = 120 # hours
dt = 0.01 # hour
Calculation = round(simulation / dt)

#  ### BOR check 하는 function 만들기
# BOG = BOG_calculation(Volume, FR, BOR, Comp_name, Pure_methane)[0]
# BOG_check = BOR_check(BOG, Comp_name, Comp_init, P_init, FR, Step, LastStep, Target, Mode_list, 30)

Mass_supply =  5000    ## kg/h
# time_range = np.arange(dt, time + dt, dt)

####  우선 압력 고정일 때,
## Mass_supply 를 dt 에 대한 배열로 입력하는것도 고려해보기

### Mode_list = ['Inerting', 'Inerting_N', 'Gassing-up', 'Cool-down', 'Cool-down_N', 'Warm-up', 'Warm-up_N', 'Aeration']
Operation = Mode_list[4]

###### Supply_condition(Comp_name, Mode_operation, Input_prop, Input1, Input2, Mass_supply, Calculation):
# Supply = Supply_condition(Comp_name, Operation, "PQmass", 0.2, 0, Mass_supply, Calculation) # CD
Supply = Supply_condition(Comp_name, Operation, "PT", 0.12, 40, Mass_supply, Calculation) # WU
# Supply = Supply_condition(Comp_name, Operation, "PT", 0.15, 20, Mass_supply, Calculation) # Inert
# Unitconv_nominal(Comp_name, Comp_init, 0.1, 40, 2, 4000) # Supply 유량 계산위한 def
# def Unitconv_nominal(Comp_name, Comp, P, T, iflag, Value):

## def Heat_transfer(Supply_condition, Init_cond, P_mode, T_ins, T_ins_ave, Ins_node, dt):

# Cooldown = Heat_transfer(Supply, Init_cond, 0, T_ins, T_ins_ave, Ins_node, dt)
# Op_inert = Heat_transfer(Supply, Init_cond, 0, T_ins, T_ins_ave, Ins_node, dt)
Warmup = Heat_transfer(Supply, Init_cond, 0, T_ins, T_ins_ave, Ins_node, dt)
# Hold = Heat_transfer(Supply, Init_cond, 1, T_ins, T_ins_ave, Ins_node, dt)

# BOR = Heat_transfer(Supply, Init_cond, 0, T_ins, T_ins_ave, Ins_node, dt)

### CHS 와 연계를 위해서 Tank info. 와 mode 같은것들은 다른걸로 빼는건 어떨지??
### 기존 계산값에서 이어서 계산할수 있는 def 도 만들기
### BOR check 하는 것 만들기...

"""""""""""""""""""""""
추출 하기
"""""""""""""""""""""""
##pd.set_option('display.float_format', '{:,.4f}'.format) 다 안먹힘
## pd.options.display.float_format = '{:.4f}'.format
## df.style.set_precision(4) # DataFrame 
Prop_col_name = ['Time', 'Tank Pressure', 'T_vap', 'T_wall', 'Density', 'Enthalpy', 'Mass_tank', 
                 'Phase', 'Qmass', 'Mass_exhaust', 'Vol_exhaust', 'Amount_spray']
Ins_col_name = np.round(Ins_node.astype(np.float64),4) # Node 소수점이 4개 고정위해서 float64로 변환
Ins_col_name = Ins_col_name.tolist()
Ins_col_name.append('Ins_ave')
Ins_col_name.append('Ins_secondary')

# df = pd.DataFrame(Cooldown[0], columns = Prop_col_name)
df = pd.DataFrame(Warmup[0], columns = Prop_col_name)

###Time 열 분리
Index_time = df['Time']

df_unit = pd.DataFrame({'Time': ['hr'],'Tank Pressure': ['MPaA'], 'T_vap': ['C'], 'T_wall': ['C'], 'Density': ['kg/m3'], 'Enthalpy': ['kJ/kg'],
   'Mass_tank': ['kg'], 'Phase': ['-'], 'Qmass': ['-'], 'Mass_exhaust': ['kg/hr'], 'Vol_exhaust': ['m3/hr'],
   'Amount_supply': ['kg']},index=['Units'])

df = pd.concat([df_unit, df])
df.set_index('Time', inplace = True)   

df_ins = pd.DataFrame(Warmup[1], columns = Ins_col_name, index = Index_time)
# df_ins.style.set_precision(4) # DataFrame
df_comp = pd.DataFrame(np.round(Warmup[2],4), columns = Comp_name_index, index = Index_time) 

"""""""""""""""""""""""
출력 하기
# """""""""""""""""""""""  
# address = 'C:\\Python_code\\04_Tank model\\Tank_result\\' 
# df.to_csv(address+'Prop.csv', header = True, index= True)          
# df_comp.to_csv(address+'Comp.csv', header = True, index= True)
# df_ins.to_csv(address+'Insulation.csv', header = True, index= True)

### GIF 만들기
# duration_rate = 0.01 # 속도
# frames_plot = []

# for plotname in plotnames:
#     if plotname.endswith(".png"):
#         # print(plotname)
#         frames_plot.append(imageio.imread(plotname))

# exportname = "Plot.gif"

# imageio.mimsave(exportname, frames_plot, format='GIF', duration=duration_rate, loop=1) # loop = 0 무한대

# for plotname in set(plotnames):
#     os.remove(plotname)