-- coding: utf-8 --

""" Spyder Editor

This is a temporary script file. """

import os, numpy as np
import math
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
Tank = T.Cargo_tank()

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

SI = RP.GETENUMdll(0,"SI WITH C").iEnum
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')    
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
            Input Data
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

## Cargo Composition (Mole base)
Comp_name =          "Nitrogen;Methane;Ethane;Propane;Butane;Isobutane;Pentane;Isopentane;Oxygen;CO2;Water"
Name_index = np.array(Comp_name.split(';'))    

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())

## Composition list
Pure_methane  = np.array([0.0,         1,    0,    0,       0,        0,     0,        0,       0,    0,   0])
Pure_nitrogen = np.array([1.0,         0,    0,    0,       0,        0,     0,        0,       0,    0,   0])
SHI_standard  = np.array([0.35,       0.88,  0.078, 0.028,  0.0105,   0,     0,        0,       0,    0,   0])
LFV_standard  = np.array([0.0047,      0.9588,  0.0304, 0.005,  0.0011,   0,     0,        0,       0,    0,   0])
BOG_standard  = np.array([0.1,         0.9,   0,    0,       0,        0,     0,        0,       0,    0,   0])
Inert_gas     = np.array([0.85  ,         0,      0,     0,       0,   0,     0,        0,    0.01,  0.14,  0])

# Comp_list = [Pure_methane, SHI_standard, LFV_standard, BOG_standard, Inert_gas]

## Atmospheric Condition
T_atm = 30 # C 외기온도
P_atm = 0.101325 # MPaA 외기압력

Voyage_day = 16    # day (pump 사양 정할때 영향)

Heat_HDC_pipe = 125 # kW (Ship loading line on deck)
# Heat_HDC_Pump = 
Heat_HDC_onshore = 7000 # kW (Onshore Pipe + pump heat)
Heat_HDC = [Heat_HDC_pipe, Heat_HDC_onshore]

"""""""""""""""""""""""""""""""""
Cargo Tank Information
"""""""""""""""""""""""""""""""""
P_tank  = 0.106             # MPaA
T_BOG   = -140                # BOG from tank

T_gassing      = 20            # Gassing-up 시 gas 온도
T_inerting     = 20            # Inerting 시 gas 온도
T_drying       = 20
T_unloading    = -140          # Unloading 시 gas 온도

N_tank = 4           # Tank 개수

## Tank 정보는 Model 표 만들어서 거기서 가져오는것으로 만들기
Circle_ratio = Tank.Circle_ratio
BOR = Tank.BOR
FR  = Tank.FR
Vol_tank_1 = Tank.Vol_tank_list[0]           # m3
Vol_tank_2 = Tank.Vol_tank_list[1]           # m3
Vol_tank_3 = Tank.Vol_tank_list[2]           # m3
Vol_tank_4 = Tank.Vol_tank_list[3]           # m3
Vol_tank_list = Tank.Vol_tank_list
Vol_tank = np.sum(Vol_tank_list)

Area_tank_1 = Tank.Area_tank_list[0]         # m2
Area_tank_2 = Tank.Area_tank_list[1]         # m2
Area_tank_3 = Tank.Area_tank_list[2]         # m2
Area_tank_4 = Tank.Area_tank_list[3]         # m2
Area_tank_list = Tank.Area_tank_list
Area_tank = np.sum(Area_tank_list)

Ins_thick = Tank.Tank_prop[0]                # mm
M_wall    = Tank.Tank_prop[1]                # kg SUS weight
Cp_wall   = Tank.Tank_prop[2]                # kJ/kg/C
M_ins     = Tank.Tank_prop[3]                # kg Ins weight
Cp_ins    = Tank.Tank_prop[4]                # kJ/kg/C
D_ins     = Tank.Tank_prop[5]                # kg/m3
Cond_ins  = Tank.Tank_prop[6]                # kJ/h/m/C
ht_NG     = Tank.Tank_prop[7]                # kJ/h/m2/C
ht_atm    = Tank.Tank_prop[8]                # kJ/h/m2/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
# BOR = 0.1           # %/day
# FR = 98.5             # %
# 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_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_list = np.array([Area_tank_1, Area_tank_2, Area_tank_3, Area_tank_4])
# Area_tank = np.sum(Area_tank_list)
# Ins_thick = 270            # mm
# M_wall = 383849            # kg SUS weight
# Cp_wall = 0.46             # kJ/kg/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

Tank_info = [Vol_tank, Area_tank, BOR, FR, Ins_thick]
Tank_prop = [M_wall, Cp_wall, M_ins, Cp_ins, D_ins, Cond_ins, ht_NG, ht_atm]
"""""""""""""""""""""""""""""""""
 Operation Criteria
"""""""""""""""""""""""""""""""""
Time_drying     = 20           # hour
Time_inerting   = 20           # hour
Time_gassing    = 20           # hour
Time_cooling    = 15           # hour
Time_loading    = 13           # hour
Time_unloading  = 12           # hour
Time_warming    = 49           # hour
Time_aerating   = 20           # hour

Operation_list = ['Time_drying', 'Time_inerting', 'Time_gassing', 'Time_cooling', 'Time_loading', 'Time_unloading', 'Time_warming', 'Time_aerating']
Time_operation = []
for i in Operation_list:
    Time_operation.append(eval(i))
# Time_operation = [Time_drying, Time_inerting, Time_gassing, Time_cooling, Time_loading, Time_unloading, Time_warming, Time_aerating]

T_CD = -130 # CD criteria
"""""""""""""""""""""""""""""""""
 Number of Equipmemt  개수와 temperature 가 함수 input 에 넣어야하는지... 이정도 input은 입력해 놔도 될지??
"""""""""""""""""""""""""""""""""
N_Cargopump = 2
N_HDC = 2

# def Pumphead (Comp_name, Comp_LNG, P_required, ):

# def Spraypump(: 

# def Heatleak():

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_tank, 0, [1.0]).Output[0]
    LH_methane = RP.REFPROPdll("Methane", "PQmass", "Heatvapz_P", SI, 1, 0, P_tank, 0, [1.0]).Output[0]

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

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

    # print(Comp_vap_LNG)
    # print(LH)
    BOG_methane = Vol * Filling_Ratio/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 MoletoMass(Name, Value): # mole to mass conversion

    Mole = []
    Comp_mass = []
    for k in Name:
        r = RP.REFPROPdll(k, 'PT', 'M', SI, 0, 0, 0.101325, 20, [1])
        Mole.append(r.Output[0])
    Sum = np.dot(Mole, Value)
    for i in range (len(Name)):
        Comp_mass.append((Mole[i]*Value[i])/Sum)

    return Comp_mass

def Comp_check(Value): ## Composition summation check

    Composition = np.array(list(map(lambda x : round(x*1000000)/1000000, Value)))
    Comp_total = np.sum(Composition)
    Comp_nor_value = Composition # Normalize 초기값
    if Comp_total == 1:
        print("Composition Complete")
    else:# normalize
        for i in range(len(Composition)):
            Comp_nor_value[i] = round((Composition[i] / Comp_total)*1000000)/1000000 

        Comp_nor_total = np.sum(Comp_nor_value)
        if Comp_nor_total == 1:
            print("Proceeding Normalize")
            print("Total : {}".format(Comp_nor_total))
        else: # Normalize 후에도 소수 4자리 이상이 되면, 가장 큰 수에 나머지를 포함시킨다.
            print("Re-Normalize")
            Max = np.max(Comp_nor_value)
            Max_index = np.where(Max == Comp_nor_value)[0][0] # [0][0] 해준 이유는 두개의 괄호안에 있으므로 위치를 찾기위해서
            Comp_nor_value[Max_index] = Comp_nor_value + (1-Comp_nor_total)
            Comp_nor_total = np.sum(Comp_nor_value)
            print("Total : {}" .format(Comp_nor_total))

    return Comp_nor_value

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 Cargopump(Time_operation, Vol_tank_list, BOG_cal, Voyage_day, N_Cargopump):

    Cargo_pump = 0
    BOG = BOG_cal[1]    # LNG 기준
    BOR = BOG_cal[2]    # LNG 기준
    Time_unloading = Time_operation[5]

    for i in Vol_tank_list:
        Arrival = i*FR/100 - (i * BOR/100 * Voyage_day)  # Vol_BOG_total = Vol_tank * BOR/100 * Voyage_day
        Capacity = Arrival / Time_unloading / N_Cargopump
        ##Unpumpable volume table 만들어야함!!!!!!!!!!!!!!!!
        ## Unloading = Arrival - unpumpable
        if Capacity > Cargo_pump :
            Cargo_pump = Capacity
        else:
            Cargo_pump = Cargo_pump

    return Cargo_pump

def HDC (Time_operation, Tank_info, Comp_name, Comp_BOG, Comp_inert, T_BOG, Overall_heat_LD, Overall_heat_CD, Prop_l, BOG_cal, Heat_HDC, N_HDC, Circle_ratio):

    Time_loading = Time_operation[5]

    Vol_tank  = Tank_info[0]
    Area_tank = Tank_info[1]
    BOR       = Tank_info[2]
    FR        = Tank_info[3]
    Ins_thick = Tank_info[4]

    BOG_LNG = BOG_cal[0]
    LH = Prop_l[0][11]       # watson's 및 HYSYS 와 차이가 많이 남.(더 작음) LH가 작아지면 BOG가 많아짐.
    P_HDC = P_tank - 0.003   # 시스템 변경 시, dP 확인필요

    #####  CASE A : Loading
    Loading_rate = Vol_tank * FR /100 / Time_loading
    # Heat_HDC_pump = Loading_rate * 
    d = RP.REFPROPdll(Comp_name, "PT", "D", SI, 1, 0, P_HDC, T_BOG, Comp_BOG).Output[0]    # BOG 밀도
    Q1 = Loading_rate * d                                                      # Q1 displaced volume of NG
    f_reduction = 30 # %, Cool-down 이후 loading 시
    Q2 = Overall_heat_LD * f_reduction /100 / LH *1000 / Time_loading             # Insulation cool-down
    Q3 = BOG_cal[1] * 0.76   # 0.76 은 50% filling 시 BOR 의 76% 적용           # BOR의 76%
    Q4 = Heat_HDC[0] * 3600 / LH                                                      # heat in leak Pipe on ship, BOG LH로 해야할지 고민중..
    Q5 = Heat_HDC[1] * 3600 / LH                                                      # Heat in leak Pipe/Pump onshore

    Q_list = np.array([Q1, Q2, Q3, Q4, Q5])
    Q_sum = np.sum(Q_list)

    HDC_loading = Q_sum / d / N_HDC
    HDC_caseA = [HDC_loading, Loading_rate, Q_list, Q_sum, d]

    #####   CASE B :Cool-down
    Required_spray_LNG = round((Overall_heat_CD*1000 / LH),-2) * 0.63 # kg
    Spray_rate = Required_spray_LNG / Time_operation[3]  # Max. Spray rate (kg/h) ==> CD 계산시 사용
    print("Spray rate :   ", Spray_rate)

    # HDC_spray = Max.(Vol_exhaust) / N_HDC   # Cool-down 을 먼저 진행해서 input 으로 받기
    # HDC_caseB = [HDC_spray, Vol_exhaust, ]

    #####  CASE C :Warm-up
    ## Refer to Warm-up Heater

    #####  CASE D :Inert gas return during gassing-up
    Exhaust_inert = Vol_tank * Circle_ratio / Time_operation[1]
    d_tank_inert = RP.REFPROPdll(Comp_name, "PT", "D", SI, 1, 0, P_tank, T_inerting, Comp_inert).Output[0]    # Inert gas 밀도
    d_HDC_inert = RP.REFPROPdll(Comp_name, "PT", "D", SI, 1, 0, P_HDC, T_inerting, Comp_inert).Output[0]    # Inert gas 밀도
    Exhaust_mass_inert = Exhaust_inert * d_tank_inert

    HDC_inert = Exhaust_mass_inert / d_HDC_inert  # m3/h, 1대로만 구동

    return HDC_caseA

def LNGvaporizer(Time_operation, Comp_name, Comp_LNG, Comp_inert, Circle_ratio, Tank_info, Cargo_pump_capacity):

    Time_gassing = Time_operation[2]
    Time_inerting = Time_operation[0]
    Tank_vol = Tank_info[0]

    ##### CASE A : Gassing-up
    d_gassing = RP.REFPROPdll(Comp_name, "PT", "D", SI, 1, 0, P_tank, T_gassing, Comp_LNG).Output[0]
    Mass_gassing = Tank_vol * Circle_ratio / Time_gassing * d_gassing  # kg/h

    h_gassing_in = RP.REFPROPdll(Comp_name, "Pliq", "h", SI, 1, 0, P_tank, 0, Comp_LNG).Output[0]
    h_gassing_out = RP.REFPROPdll(Comp_name, "PT", "h", SI, 1, 0, P_tank, T_gassing, Comp_LNG).Output[0]
    Power_gassing = Mass_gassing * (h_gassing_out - h_gassing_in) / 3600  # kW

    LNGV_caseA = [Mass_gassing, Power_gassing]

    ##### CASE B : Unloading
    CP_capa = Cargo_pump_capacity
    Total_N_CP = N_Cargopump * N_tank                                         # C/P 총 개수
    Vol_max = CP_capa * Total_N_CP                                             # 최대 unloading volume flowrate

    d_unloading = RP.REFPROPdll(Comp_name, "PT", "D", SI, 1, 0, P_tank, T_unloading, Comp_LNG).Output[0]
    Mass_unloading = Vol_max * d_unloading

    h_unloading_in = RP.REFPROPdll(Comp_name, "Pliq", "h", SI, 1, 0, P_tank, 0, Comp_LNG).Output[0]
    h_unloading_out = RP.REFPROPdll(Comp_name, "PT", "h", SI, 1, 0, P_tank, T_unloading, Comp_LNG).Output[0]
    Power_unloading = Mass_unloading * (h_unloading_out - h_unloading_in) / 3600  # kW

    LNGV_caseB = [Mass_unloading, Power_unloading]

    ##### CASE C : Inerting with LN2
    d_inerting = RP.REFPROPdll(Comp_name, "PT", "D", SI, 1, 0, P_tank, T_inerting, Comp_nitrogen).Output[0]
    print(d_inerting)
    Mass_inerting = Tank_vol * Circle_ratio / Time_inerting * d_inerting  # kg/h

    h_inerting_in = RP.REFPROPdll(Comp_name, "Pliq", "h", SI, 1, 0, P_tank, 0, Comp_nitrogen).Output[0]
    h_inerting_out = RP.REFPROPdll(Comp_name, "PT", "h", SI, 1, 0, P_tank, T_inerting, Comp_nitrogen).Output[0]
    Power_inerting = Mass_inerting * (h_inerting_out - h_inerting_in) / 3600  # kW

    LNGV_caseC = [Mass_inerting, Power_inerting]
    return LNGV_caseA, LNGV_caseB, LNGV_caseC

def IGG(Time_operation, Comp_name, Comp_inert, Tank_info, Circle_ratio):

    Vol_tank = Tank_info[0]
    Time_drying = Time_operation[0]
    Time_inerting = Time_operation[1]
    Time_aeration = Time_operation[7]
    T_operation = 15     # 왜 15도 인지는 잘 모르겠음..

    ##### CASE A : DRYING
    d_drying = RP.REFPROPdll("Air", "PT", "D", SI, 1, 0, P_atm, T_operation,[1.0] ).Output[0]
    d_drying_nominal = RP.REFPROPdll("Air", "PT", "D", SI, 1, 0, P_atm, 0 ,[1.0] ).Output[0]

    IGG_drying = Vol_tank * Circle_ratio / Time_drying * d_drying / d_drying_nominal

    ##### CASE B : INERTING BEFORE GAS FILLING
    d_inerting = RP.REFPROPdll(Comp_name, "PT", "D", SI, 1, 0, P_atm, T_operation, Comp_inert ).Output[0]
    d_inerting_nominal = RP.REFPROPdll(Comp_name, "PT", "D", SI, 1, 0, P_atm, 0 , Comp_inert ).Output[0]

    IGG_inerting = Vol_tank * Circle_ratio / Time_inerting * d_inerting / d_inerting_nominal    

    ##### CASE C : INERTING AFTER WARMING-UP
    T_after_WU = 5
    d_inerting_after_WU = RP.REFPROPdll(Comp_name, "PT", "D", SI, 1, 0, P_atm, T_after_WU, Comp_inert ).Output[0]

    IGG_inerting_after_WU = Vol_tank * Circle_ratio / Time_inerting * d_inerting_after_WU / d_inerting_nominal

    ##### CASE D : AERATION
    d_aerating = RP.REFPROPdll("Air", "PT", "D", SI, 1, 0, P_atm, T_operation,[1.0] ).Output[0]
    d_aerating_nominal = RP.REFPROPdll("Air", "PT", "D", SI, 1, 0, P_atm, 0 ,[1.0] ).Output[0]

    IGG_aerating = Vol_tank * Circle_ratio / Time_aerating * d_aerating / d_aerating_nominal        

    return IGG_drying, IGG_inerting, IGG_inerting_after_WU 

"""""""""""""""""""""""""""""""""
   Calculation
"""""""""""""""""""""""""""""""""

Comp_mole = LFV_standard                                                       # Composition setting
Comp_mass = MoletoMass(Name_index,Comp_mole)                                   # Mole to Mass conversion
Comp_LNG = Comp_check(Comp_mass)                                               # Composition sum 1 check

Comp_mole = BOG_standard                                                       # Composition setting
Comp_mass = MoletoMass(Name_index,Comp_mole)                                    # Mole to Mass conversion
Comp_BOG = Comp_check(Comp_mass)                                               # Composition sum 1 check

Comp_mole = Inert_gas
Comp_mass = MoletoMass(Name_index,Comp_mole)
Comp_inert = Comp_check(Comp_mass)

Comp_nitrogen = Pure_nitrogen
### Composition 결정 후 계산되는 값
Prop_l = Properties(Comp_name, Comp_LNG, "PQmass", P_tank, 0)

T_boiling = RP.REFPROPdll(Comp_name, "Pliq", "T", SI, 1, 0, P_tank, 0, Comp_LNG).Output[0]
Overall_heat_LD = ( M_wall * Cp_wall * (T_atm - T_boiling) + M_ins * Cp_ins * (T_atm - (T_atm + T_boiling)/2) )/1000 # MJ
Overall_heat_CD = ( M_wall * Cp_wall * (T_atm - T_CD) + M_ins * Cp_ins * (T_atm - (T_atm + T_CD)/2) )/1000 # MJ
print("Overall heat Capacity : {}" .format(Overall_heat_LD))
print("Overall heat for Cool-down : {}" .format(round(Overall_heat_CD*1000 / Prop_l[0][11],2) * 0.63 / Time_operation[3]))

BOG_cal = BOG_calculation(Vol_tank, FR, BOR, Comp_name, Comp_LNG)

Cargo_pump_capacity = Cargopump(Time_operation, Vol_tank_list, BOG_cal, Voyage_day, N_Cargopump)
HDC_capacity = HDC (Time_operation, Tank_info, Comp_name, Comp_BOG, Comp_inert, T_BOG, Overall_heat_LD, Overall_heat_CD, Prop_l, BOG_cal, Heat_HDC, N_HDC, Circle_ratio)
LNGV_capacity = LNGvaporizer(Time_operation, Comp_name, Comp_LNG, Comp_inert, Circle_ratio, Tank_info, Cargo_pump_capacity)
IGG_capacity = IGG(Time_operation, Comp_name, Comp_inert, Tank_info, Circle_ratio)

def Properties(Comp_name, Comp_value, Input_prop, input1, input2): Prop_list = "P;T;D;H;M;CP;Phase;Qmass;VIS;TCX;PRANDTL;KV;Heatvapz_P" # VISCOSITY 는 100만 나누어야 함. Prop_index = np.array(Prop_list.replace(';',' ').split()) Prop = [] Value = [] for k in Prop_index: # P;T;D;H;M;CP;VIS;TCX;PRANDTL;KV;Heatvapz_P r = RP.REFPROPdll(Comp_name, Input_prop , k, SI, 0,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

`# -*- coding: utf-8 -*- """ Spyder Editor This is a temporary script file. """ """ 22. 03. 17 """ import os, numpy as np import math 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 import scipy.integrate as int 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;Phase;Qmass;VIS;TCX;PRANDTL;KV;Heatvapz_P" # VISCOSITY 는 100만 나누어야 함. Prop_index = np.array(Prop_list.replace(';',' ').split()) # df1 = pd.read_pickle('Pipetable.p') ''''''''''''''''''''''''''' Initial Condition ''''''''''''''''''''''''''' Comp_name = "Nitrogen;Methane;Ethane;Propane;Butane;Isobutane;Pentane;Isopentane;Oxygen;CO2;Water" Comp_init = np.array([0.0, 100.0, 0, 0, 0, 0, 0, 0, 0, 0, 0]) T_atm = 45 # C P_init = 0.11 # Mpa Volume = 14012 # m3 Surface = 3365 # m2 Filling = 98 # % Vol_liq = Volume * Filling/100 Vol_vap = Volume - Vol_liq ''''''''''''''''''''''''''' Equilibrium Mode ''''''''''''''''''''''''''' Qmass = 0.00010902809999999928 # initial value Step = 0.01 LastStep = 0.000000000001 Target = 0 direction = 0 D_tank = RP.REFPROPdll(Comp_name,"PQmass","D",SI,0,0,P_init,Qmass,Comp_init).Output[0] D_liq = RP.REFPROPdll(Comp_name,"PQmass","Dliq",SI,0,0,P_init,Qmass,Comp_init).Output[0] D_vap = RP.REFPROPdll(Comp_name,"PQmass","Dvap",SI,0,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 # 질량보존 확인항 if Mass_conv < Target: direction = 1 else: direction = -1 if Mass_conv == Target: sys.exit(0) # sys.exit(0) 프로그램 정상 종료, (1) 강제 종료 if Filling == 0: # Qmass = 1 T_vap = T_atm T_lng = T_atm T_liq = T_atm T_wall = T_atm T_ins = T_atm # Filling 이 0일때 property, compostion 계산하는 항 만들기. else: for i in range(1,10000,1): Qmass = Qmass + Step * direction if Mass_conv < Target: if direction == 1: if Step > LastStep: Step = Step / 10 direction = direction * -1 # else: # exit() else: if direction == -1: if Step > LastStep: Step = Step / 10 direction = direction * -1 # else: # exit() if Step < LastStep: break D_tank = RP.REFPROPdll(Comp_name,"PQmass","D",SI,0,0,P_init,Qmass,Comp_init).Output[0] D_liq = RP.REFPROPdll(Comp_name,"PQmass","Dliq",SI,0,0,P_init,Qmass,Comp_init).Output[0] D_vap = RP.REFPROPdll(Comp_name,"PQmass","Dvap",SI,0,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 # pri``nt (Mass_tank, Mass_liq, Mass_vap, D_liq, D_tank, D_vap) Mass_conv = Mass_tank-Mass_liq-Mass_vap # print (Qmass) Comp_liq_mole = np.array(RP.REFPROPdll(Comp_name,"PQmass","xmoleliq",SI,0,0,P_init,Qmass,Comp_init).Output[0:Comp_name.count(';')]) Comp_vap_mole = np.array(RP.REFPROPdll(Comp_name,"PQmass","xmolevap",SI,0,0,P_init,Qmass,Comp_init).Output[0:Comp_name.count(';')]) Pro_lng_mole = np.array(RP.REFPROPdll(Comp_name,"PQmass",Prop_list,SI,0,0,P_init,0,Comp_init).Output[0:Prop_list.count(';')]) Pro_liq_mole = np.array(RP.REFPROPdll(Comp_name,"Pliq",Prop_list,SI,0,0,P_init,0,Comp_liq_mole).Output[0:Prop_list.count(';')]) Pro_vap_mole = np.array(RP.REFPROPdll(Comp_name,"Pvap",Prop_list,SI,0,0,P_init,0,Comp_vap_mole).Output[0:Prop_list.count(';')]) T_vap = Pro_vap_mole[0] T_wall = T_vap # 초기 조건 ht_NG = 20.93 # kJ/h/m2/C ht_atm = 10.47 # kJ/h/m2/C ''''''''''''''''''''''''''' Select Mode ''''''''''''''''''''''''''' Mode_list = ['Inerting', 'Inerting_N', 'Gassing-up', 'Cool-down', 'Cool-down_N', 'Warm-up', 'Warm-up_N', 'Aeration'] Mode_pressure = "Constant" # Constant, Variable ''''''''''''''''''''''''''' Insulation information ''''''''''''''''''''''''''' M_wall = 383849 # kg Cp_wall = 0.46 # kJ/kg/C T_wall = Pro_lng_mole[0] # C M_ins = 1891303 # kg Cp_ins = 1.07 # kJ/kg/C D_ins = 190 # kg/m3 Cond_ins = 0.0963 # kJ/h/m/C Ins_thick = 270 # mm ''''''''''''''''''''''''''''''''''''''' Insulation Initial condition ''''''''''''''''''''''''''''''''''''''' n = 12 # Insulation node 수 len_interval = Ins_thick/1000 / n # Node 별 거리 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/len_interval) # ins del_T # step_node = np.arange(start = 0, stop = Ins_thick/1000 + 0.0000001, step = Ins_thick/n/1000) T_ins = np.zeros((1,n+1)) # Insulation 초기온도 조건 if T_ins_step == 0: T_ins = np.full((1,n+1), T_ins) else: T_ins = np.linspace(T_wall, T2, n+1) T_ins = T_ins.T # df_Tins = pd.DataFrame(T_ins, columns =step_node) # Column을 step_node 로 #df_ins = df_ins.rename(columns=df_ins.iloc[0]) # 첫 행을 columns name 으로 dt_ins = 0 diffusivity = Cond_ins / (Cp_ins * D_ins) # 제공 초기조건 Input_prop = "PQmass" # 2가지 고르기 (input에서) Input1 = 0.2 # MPaA Input 1 값 Input2 = 0 # Input 2 값 # Supply 조건 지정. # ### BOR check 하는 function 만들기 time = 12 # hours t_step = 0.01 # hour def Properties(Comp_name, Comp_value, Input_prop, input1, input2): Prop = [] Value = [] for k in Prop_index: # P;T;D;H;M;CP;VIS;TCX;PRANDTL;KV;Heatvapz_P r = RP.REFPROPdll(Comp_name, Input_prop , k, SI, 0,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 # Mode_list = ['Inerting', 'Inerting_N', 'Gassing-up', 'Cool-down', 'Cool-down_N', 'Warm-up', 'Warm-up_N', 'Aeration'] def Supply_condition(Comp_name, Mode_operation, Input_prop, Input1, Input2): if Mode_operation == Mode_list[0]: Comp_supply = np.array([85.0, 0, 0, 0, 0, 0, 0, 0, 1, 14, 0]) elif Mode_operation == Mode_list[1] or Mode_operation == Mode_list[4] or Mode_operation == Mode_list[6]: Comp_supply = np.array([100.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]: Comp_supply = Comp_init else: # Aeration Comp_supply = np.array([78.0, 0, 0, 0, 0, 0, 0, 0, 21, 1, 0]) Value = Properties(Comp_name, Comp_supply, Input_prop, Input1, Input2)[0] return Comp_supply, Value #### 우선 압력 고정일 때, ## Mass_supply 를 dt 에 대한 배열로 입력하는것도 고려해보기 def Heat_transfer(Comp_tank, Comp_supply, Mass_supply, P_mode): dh = (-1 * (Mass_supply * (h - h_supply) + ht_NG * Surface * (T_vap - T_wall)) * dt) / Mass_tank dT_wall = ((ht_NG * Surface * (T_vap - T_wall) + ht_atm * Surface_area * (T_atm - T2)) * dt - M_ins * Cp_ins* T_ins_var) / (M_wall * Cp_wall) h = h + dh T_wall = T_wall + dT_wall Mass_exhaust_init = 100 # 100 kg/h # 아래 항을 mole 로 바로 계산할 수 있으면 좋을텐데? Comp_init 은 mass base # Comp = (Comp_init * Tank_mass + Comp_supply * Mass_supply * dt) / (Tank_mass) Comp = list(map(lambda x, y: (x * Mass_tank + y * Mass_supply * dt) / (Mass_tank + Mass_supply*dt), Comp, Comp_supply)) # 압력항에 대한 것도 같이 넣어줘야 할것 같음. if P_mode == 0: # 0: Const, 1: Various d = RP.REFPROPdll(Comp_name, "PH" , "D", SI, 0,0, P_init, h, Comp).Output[0] else: d = D_tank # Mass_conv = Mass_supply * dt + Mass_tank - Mass_exhaust Mass_tank = Volume * d # for i in len() # ins_T[i] = Thermal diffusivity x dt / Interval length * (ins 거리0 - ins 거리2) + (1- 2 x (Thermal diffusivity x dt / Interval length **2)) x Ins 거리 1 # for i in range (0,12/3600+t_step,t_step): Supply = Supply_condition(Comp_name, Mode_list[1], Input_prop, Input1, Input2) Prop = Properties(Comp_name, Comp_init, "PT", 0.12, 0) `

# -*- coding: utf-8 -*- """ Spyder Editor This is a temporary script file. """ """ 22. 06. 27 """ import os, numpy as np import math 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 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') ''''''''''''''''''''''''''' Initial Condition ''''''''''''''''''''''''''' Comp_name = "Nitrogen;Methane;Ethane;Propane;Butane;Isobutane;Pentane;Isopentane;Oxygen;CO2;Water" Comp_init = 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 T_atm = 40 # C P_init = 0.13 # Mpa T_set_vapor = 40 Volume = 18000 # m3 Surface = 4200 # m2 Filling = 0 # % Vol_liq = Volume * Filling/100 Vol_vap = Volume - Vol_liq ''''''''''''''''''''''''''''' Insulation information ''''''''''''''''''''''''''''' M_wall = 60205 # kg Cp_wall = 0.35 # kJ/kg/C # T_wall = T_vap # 초기 조건 M_ins = 424071 # kg Cp_ins = 1.07 # kJ/kg/C D_ins = 118 # kg/m3 Cond_ins = 0.08 # kJ/h/m/C Ins_thick = 400 # mm ht_NG = 37.46 # kJ/h/m2/C ht_atm = 16.8 # kJ/h/m2/C ''''''''''''''''''''''''''' Equilibrium Mode ''''''''''''''''''''''''''' # Qmass = 0.1 # initial value funcion 에 기입 Step = 0.01 LastStep = 0.000000000001 Target = 0.001 # direction = 0 ''''''''''''''''''''''''''' Plot ''''''''''''''''''''''''''' x_axis = [] y_axis = [] y_axis2 = [] y_axis3 = [] ############################################################## ####################### Definition ########################### ############################################################## # Mode_list = ['Inerting', 'Inerting_N', 'Gassing-up', 'Cool-down', 'Cool-down_N', 'Warm-up', 'Warm-up_N', 'Aeration'] # Input_prop = 'PT', 'PH' 등등 def Supply_condition(Comp_name, Mode_operation, Input_prop, Input1, Input2, Mass_supply, Calculation): 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]: 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, Surface, Filling, Vol_liq, Vol_vap, Step, LastStep, Target): Qmass = 0.1 # 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] # (hFld, hIn, hOut, iUnits, iMass, iFlag, a, b, z, Output, hUnits, iUCode, x, y, x3, q, ierr, herr, hFld_length, hIn_length, hOut_length, hUnits_length, herr_length) 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 Mass_conv < Target: ### Qmass 가 음수가 되는 경우 피하기 위한 함수 if Qmass < 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 # print(Qmass, Mass_conv, D_tank, D_liq, D_vap) 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 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 T_vap = Init_condition(Comp_name, Comp_init, P_init, Surface, Filling, Vol_liq, Vol_vap, Step, LastStep, Target)[4][1] T_lng = Init_condition(Comp_name, Comp_init, P_init, Surface, Filling, Vol_liq, Vol_vap, Step, LastStep, Target)[2][1] T_liq = Init_condition(Comp_name, Comp_init, P_init, Surface, Filling, Vol_liq, Vol_vap, Step, LastStep, Target)[3][1] T_wall = T_vap Mass_tank = Volume * Init_condition(Comp_name, Comp_init, P_init, Surface, Filling, Vol_liq, Vol_vap, Step, LastStep, Target)[2][2] Mass_liq = Volume * Init_condition(Comp_name, Comp_init, P_init, Surface, Filling, Vol_liq, Vol_vap, Step, LastStep, Target)[3][2] Mass_vap = Volume * Init_condition(Comp_name, Comp_init, P_init, Surface, Filling, Vol_liq, Vol_vap, Step, LastStep, Target)[4][2] ''''''''''''''''''''''''''' Select Mode ''''''''''''''''''''''''''' Mode_list = ['Inerting', 'Inerting_N', 'Gassing-up', 'Cool-down', 'Cool-down_N', 'Warm-up', 'Warm-up_N', 'Aeration'] 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,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 # step_node = np.arange(start = 0, stop = Ins_thick/1000 + 0.0000001, step = Ins_thick/n/1000) T_ins = np.zeros(n+1) # Insulation 초기온도 조건 if T_ins_step == 0: T_ins = np.full(n+1, T_wall) else: T_ins = np.linspace(T_wall, T2, n+1) # T_ins = T_ins.T 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) # df_Tins = pd.DataFrame(T_ins, columns =step_node) # Column을 step_node 로 #df_ins = df_ins.rename(columns=df_ins.iloc[0]) # 첫 행을 columns name 으로 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 def Heat_transfer(Supply_condition, Init_cond, P_mode, T_ins, T_ins_ave, Ins_node): #################### 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[4][1] h = Init_cond[4][3] d = Init_cond[4][2] Phase = Init_cond[5] Qmass = Init_cond[2][6] Comp = Comp_init Comp_liq = Init_cond[0] Comp_vap = Init_cond[1] Supply = Supply_condition T_wall = T_vap Mass_tank = Volume * Init_cond[2][2] T_ins_var = 0 T_ins = T_ins Ins_node = Ins_node T_ins_temp = np.zeros(len(T_ins)) ################### Data save ############################################# T_ins_table = np.append(T_ins,T_ins_ave) # print(T_ins_table) Prop_table = np.array([P_init, T_vap, T_wall, d, h, Mass_tank, Phase, Qmass]) Comp_table = np.array([Comp_init, Comp_liq, Comp_vap]) # print(Comp_table) for i in range (Calculation): # if Operation == "Cool-down": # if T_vap < -140: # break # elif Operation == "Cool-down": dh = (-1 * (Supply[2][i] * (h - Supply[1][3]) + ht_NG * Surface * (T_vap - T_wall)) * dt) / Mass_tank dT_wall = ((ht_NG * Surface * (T_vap - T_wall) + ht_atm * Surface * (T_atm - T2)) * dt - M_ins * Cp_ins* T_ins_var) / (M_wall * Cp_wall) # h = float(int((h + dh)*1000)/1000) h = h + dh T_wall = T_wall + dT_wall Mass_exhaust_init = 100 # 100 kg/h # 아래 항을 mole 로 바로 계산할 수 있으면 좋을텐데? Comp_init 은 mass base # Comp = (Comp_init * Tank_mass + Comp_supply * Mass_supply * dt) / (Tank_mass) Comp = list(map(lambda x, y: (x * Mass_tank + y * Supply[2][i] * dt) / (Mass_tank + Supply[2][i]*dt), Comp, Supply[0])) # Comp = np.array(Comp) # 압력항에 대한 것도 같이 넣어줘야 할것 같음. if P_mode == 0: # 0: Const, 1: Various d = RP.REFPROPdll(Comp_name, "PH" , "D", SI, 1, 0, P_init, h, Comp).Output[0] Mass_tank = Volume * d P = P_init else: Mass_tank = Mass_tank + Supply[2][i]*dt d = Mass_tank / Volume P = RP.REFPROPdll(Comp_name, "DH" , "P", SI, 1, 0, d, h, Comp).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(';')]) T_vap = np.array(RP.REFPROPdll(Comp_name,"DH","T",SI,1,0,d,h,Comp).Output[0]) Qmass = np.array(RP.REFPROPdll(Comp_name,"DH","Qmass",SI,1,0,d,h,Comp).Output[0]) Phase = np.array(RP.REFPROPdll(Comp_name,"DH","Phase",SI,1,0,d,h,Comp).hUnits) print(i*dt, dh, T_vap, T_wall) Comp_liq_mass = np.array(RP.REFPROPdll(Comp_name,"DH","xmassliq",SI,1,0,P_init,Qmass,Comp).Output[0:Comp_name.count(';')]) Comp_vap_mass = np.array(RP.REFPROPdll(Comp_name,"DH","xmassvap",SI,1,0,P_init,Qmass,Comp).Output[0:Comp_name.count(';')]) Comp_table = np.vstack([Comp_table, [Comp, Comp_liq, Comp_vap]]) # Prop_table = np.array([P, Prop[1], d, h, T_wall, Mass_tank, 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_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_table = np.vstack([T_ins_table, T_ins_temp_table]) # Time step 별 Insul node 별 테이블화 # Prop_table = np.array([P_init, T_vap, T_wall, d, h, Mass_tank, Phase, Qmass]) Prop_tank = np.array([P, T_vap, T_wall, d, h, Mass_tank, Phase, Qmass]) 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 ################################ 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) 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) 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): # Prop = Properties(Comp_name, Comp_init, "PT", 0.12, 0) """"""""""""""""""""""" Calculation """"""""""""""""""""""" # Supply 조건 지정. Input_prop = "PQmass" # 2가지 고르기 (input에서) Input1 = 0.2 # MPaA Input 1 값 Input2 = 0 # Input 2 값 # ### BOR check 하는 function 만들기 time = 17 # hours dt = 0.05 # hour Calculation = round(time / dt) Mass_supply = 2800 ## 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[3] Supply = Supply_condition(Comp_name, Operation, "PQmass", 0.128, 0, Mass_supply, Calculation) Init_cond = Init_condition(Comp_name, Comp_init, P_init, Surface, Filling, Vol_liq, Vol_vap, Step, LastStep, Target) # Heat_transfer(Supply_condition, Init_cond, Mass_supply, P_mode, T_ins, T_ins_ave, Ins_node): ## 계산 Heat_transfer(Supply, Init_cond, 0, T_ins, T_ins_ave, Ins_node) Cooldown = Heat_transfer(Supply, Init_cond, 0, T_ins, T_ins_ave, Ins_node)

JungHulk / Hulk-Engineering

eq des #1

-- coding: utf-8 --