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.

Sign up for GitHub

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

Jump to bottom

eq des #1

Open
JungHulk opened this issue Mar 17, 2022 · 5 comments
Open

eq des #1

JungHulk opened this issue Mar 17, 2022 · 5 comments

Comments

@JungHulk
Copy link
Owner

JungHulk commented Mar 17, 2022 •
edited
Loading

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

https://trc.nist.gov/refprop/10-BETA/REFPROP.HTM

@JungHulk

This comment was marked as off-topic.

Sign in to view
@JungHulk
Copy link
Owner Author

JungHulk commented Mar 18, 2022 •
edited
Loading 

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

@JungHulk
Copy link
Owner Author

JungHulk commented Mar 25, 2022 •
edited
Loading 

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

@JungHulk
Copy link
Owner Author

JungHulk commented Jul 6, 2022 •
edited
Loading 

# -*- 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 JungHulk changed the title Tank modeling-1 eq des Sep 7, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
1 participant
@JungHulk