Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Basic HTE #11

Open
JungHulk opened this issue Feb 1, 2024 · 1 comment
Open

Basic HTE #11

JungHulk opened this issue Feb 1, 2024 · 1 comment

Comments

@JungHulk
Copy link
Owner

JungHulk commented Feb 1, 2024

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
Copy link
Owner Author

JungHulk commented Feb 1, 2024

# -*- 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)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant