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

Tank_operation #3

Open
JungHulk opened this issue Aug 1, 2022 · 2 comments
Open

Tank_operation #3

JungHulk opened this issue Aug 1, 2022 · 2 comments

Comments

@JungHulk
Copy link
Owner

JungHulk commented Aug 1, 2022 •
edited
Loading

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

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])
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_atmIns_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)

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

JungHulk commented Sep 7, 2022

information of vessel.

Calling the Class for operation.

-- 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 Sep 7, 2022

-- coding: utf-8 --

"""
Created on Wed Aug 10 13:36:23 2022

@author: 252914
"""

-- coding: utf-8 --

"""
22. 08. 10.

  1. BOR base kW 맞도록 insulation conductivity 조절
  2. BOR base massflowrate 맞도록 LNG ocnvective heat transfer 조절
  3. Holding time 계산

"""

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

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, T_wall):

#################### 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%로 진행
dt = dt
# dtt = 2 * dt
T_wall = T_wall
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
P = RP.REFPROPdll(Comp_name, "DH" , "P", SI, 1, 1, d, h, Comp).Output[0]
# print("=============P========= : ", P)


################### 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])
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
        # h= (Mass_liq * h_liq + (Mass_vap - Mass_exhaust * dt) * h_vap) / Mass_tank
        # print(((ht_atm * Surface) * (T_atm - T_ins[-1]))/509)
        # print(((ht_NG * Surface * (T_vap - T_wall)) * dt), Mass_liq * h_liq, Mass_vap*h_vap)


        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
    fill = Mass_liq / d_liq / Volume * 100
    
    print('Time :', "%-8s" % np.round(i*dt+dt,2), 'Press. :', "%-8s" % np.round(P,6), 'T_vap :', "%-8s" % np.round(T_vap,3), 'T_wall :', "%-8s" % np.round(T_wall,2), 'T_sec :', np.round(T_ins_sec,2) )
    print('Filling Ratio :  ', np.round(fill,4))
    
    Heat_to_gas = Mass_tank * (h - h_pre) / dt        # kJ/h
    Heat_to_pri = M_wall * Cp_wall * T_ins_var / dt   # kJ/h
    Heat_total  = U * Surface * (T_atm - T_vap)     # kJ/h
    BOG = Heat_total / Supply[1][11] # kg/h
    
    # print(Heat_to_gas, Heat_to_pri, Heat_total)
    if i == 1/dt:
        print("Heat_total", Heat_total,"Cal_BOG :", BOG)
        time.sleep(1)
        break
        
    
    
    ################################ 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(-160,-155)
    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]

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

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,1000,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,1,P_init,Qmass,Comp_init).Output[0:Comp_name.count(';')])
    Comp_vap_mass = np.array(RP.REFPROPdll(Comp_name,"PQmass","xmassvap",SI,1,1,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
    
return Nominal, Vol, Mass

def BOG_calculation(Tank_volume, Filling_Ratio, BOR, Comp_name, Comp_LNG, U, Area, T_vap, T_atm): # 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

Heat_leak = (U * Surface * (T_atm-T_vap))/3600

print("BOG Methane :", "%-5s" % round(BOG_methane,2), "%-5s" % "kg/h")
print("BOG_LNG :", "%-5s" % round(BOG_LNG,2), "kg/h")
print("Heat leak :", "%-5s" % round(Heat_leak,2), "kW")
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([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
T_atm = 45 # C
P_init = 0.15 # Mpa
T_set_vapor = -150

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

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 각 위치
for i in range (1000):
print(i)
T2 = (T_wall + (ht_atm * Ins_thick/1000/Cond_ins)T_atm) / (1+(ht_atmIns_thick /1000/Cond_ins))
T_wall_update = (T2 + (ht_NG * Ins_thick/1000/Cond_ins)T_vap) / (1+(ht_NGIns_thick /1000/Cond_ins))
check = T_wall - T_wall_update
print(check, T2, T_wall, T_wall_update)
if abs(check) < 0.1:
break
else:
T_wall = T_wall_update

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

"""""""""""""""""""""""
Calculation
"""""""""""""""""""""""
simulation = 1*24 # hours
dt = 0.2 # hour
Calculation = round(simulation / dt)
Mass_supply = 0 ## kg/h
Operation = Mode_list[2]

Supply_condition(Comp_name, Mode_operation, Input_prop, Input1, Input2, Mass_supply, Calculation):

Zero = Supply_condition(Comp_name, Operation, "PQmass", 0.13, 0, Mass_supply, Calculation)
Hold = Heat_transfer(Zero, Init_cond, 1, T_ins, T_ins_ave, Ins_node, dt, T_wall)

Check = BOG_calculation(Volume, FR, BOR, Comp_name, Comp_init, U, Surface, T_vap, T_atm)

return BOG_methane, BOG_LNG, Revised_BOR

Hold = Heat_transfer(Zero, Init_cond, 1, T_ins, T_ins_ave, Ins_node, dt, T_wall)

"""""""""""""""""""""""

추출 하기

"""""""""""""""""""""""

##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_supply']

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(Hold[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(Hold[1], columns = Ins_col_name, index = Index_time)

df_ins.style.set_precision(4) # DataFrame

df_comp = pd.DataFrame(np.round(Hold[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
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
1 participant
@JungHulk