search

JungHulk / Hulk-Engineering

1 stars 0 forks source link

Tank_operation #3

Open JungHulk opened 2 years ago

JungHulk commented 2 years ago

-- coding: utf-8 --

"""

    1. 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_atm*Ins_thick /1000/Cond_ins))

Ins Outer temp.

T_ins_step = (T2 - T_wall) / (Ins_thick/1000/Ins_interval) # ins del_T T_ins = np.zeros(n+1)

Insulation secondary barrier 위치 정하기

Ins_sec = Tank.sec if np.any(Ins_node == Ins_sec) == True: # Insulation node 와 secondary barrier 가 같은경우, loc_sec = np.where(Ins_node == Ins_sec)[0][0] #Insulation location loc_sec_1 = np.where(Ins_node == Ins_sec)[0][0]

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

else: # 다를 경우는, 근처 2개 값의 평균값을 사용 nearest = Ins_node.flat[np.abs(Ins_node - Ins_sec).argmin()] if nearest > Ins_sec: Ins_sec_1 = Ins_sec - 0.01 elif nearest < Tank.sec: Ins_sec_1 = Ins_sec + 0.01 nearest_1 = Ins_node.flat[np.abs(Ins_node - Ins_sec_1).argmin()] loc_sec = np.where(Ins_node == nearest)[0][0] loc_sec_1 = np.where(Ins_node == nearest_1)[0][0]

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

Insulation 초기온도 조건

if T_ins_step == 0: T_ins = np.full(n+1, T_wall) else: T_ins = np.linspace(T_wall, T2, n+1)

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

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

return integration

T_ins_ave = simpson(T_ins,n)

dt_ins = 0

diffusivity = Cond_ins / (Cp_ins * D_ins)

제어기? 혹은 기준설정 후 자동으로 종료되게 만들기

def Operation_criteria(Mode_list, Operation, Temperature, Composition):

[0] Inerting

[1] Inerting_N

[2] Gassing-up

[3] Cool-down

[4] Cool-down_N

[5] Warm-up

[6] Warm-up_N

[7] Aeration

Operation = Mode_list

if Operation.find("Cool") is not -1:

if Temperature < -140:

break

elif Operation.find("Inert") is not -1:

if Composition < 0.02:

break

Prop = Properties(Comp_name, Comp_init, "PT", 0.12, 0)

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

Supply 조건 지정.

Input_prop = "PQmass" # 2가지 고르기 (input에서)

Input1 = 0.106 # MPaA Input 1 값

Input2 = 0 # Input 2 값

simulation = 120 # hours dt = 0.01 # hour Calculation = round(simulation / dt)

BOR check 하는 function 만들기

BOG = BOG_calculation(Volume, FR, BOR, Comp_name, Pure_methane)[0]

BOG_check = BOR_check(BOG, Comp_name, Comp_init, P_init, FR, Step, LastStep, Target, Mode_list, 30)

Mass_supply = 5000 ## kg/h

time_range = np.arange(dt, time + dt, dt)

우선 압력 고정일 때,

Mass_supply 를 dt 에 대한 배열로 입력하는것도 고려해보기

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

Operation = Mode_list[4]

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

Supply = Supply_condition(Comp_name, Operation, "PQmass", 0.2, 0, Mass_supply, Calculation) # CD

Supply = Supply_condition(Comp_name, Operation, "PT", 0.12, 40, Mass_supply, Calculation) # WU

Supply = Supply_condition(Comp_name, Operation, "PT", 0.15, 20, Mass_supply, Calculation) # Inert

Unitconv_nominal(Comp_name, Comp_init, 0.1, 40, 2, 4000) # Supply 유량 계산위한 def

def Unitconv_nominal(Comp_name, Comp, P, T, iflag, Value):

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

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

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

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

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

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

CHS 와 연계를 위해서 Tank info. 와 mode 같은것들은 다른걸로 빼는건 어떨지??

기존 계산값에서 이어서 계산할수 있는 def 도 만들기

BOR check 하는 것 만들기...

""""""""""""""""""""""" 추출 하기 """""""""""""""""""""""

pd.set_option('display.float_format', '{:,.4f}'.format) 다 안먹힘

pd.options.display.float_format = '{:.4f}'.format

df.style.set_precision(4) # DataFrame

Prop_col_name = ['Time', 'Tank Pressure', 'T_vap', 'T_wall', 'Density', 'Enthalpy', 'Mass_tank', 'Phase', 'Qmass', 'Mass_exhaust', 'Vol_exhaust', 'Amount_spray'] Ins_col_name = np.round(Ins_node.astype(np.float64),4) # Node 소수점이 4개 고정위해서 float64로 변환 Ins_col_name = Ins_col_name.tolist() Ins_col_name.append('Ins_ave') Ins_col_name.append('Ins_secondary')

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

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 commented 2 years ago

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 commented 2 years ago

-- coding: utf-8 --

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

@author: 252914 """

-- coding: utf-8 --

"""

    1. 10.

  2. BOR base kW 맞도록 insulation conductivity 조절

  3. BOR base massflowrate 맞도록 LNG ocnvective heat transfer 조절

  4. 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)