search

JungHulk / Hulk-Engineering

1 stars 0 forks source link

SJ_PRV #4

Open JungHulk opened 2 years ago

JungHulk commented 2 years ago

import os, numpy as np import math import matplotlib.pyplot as plt

from win32com.client import Dispatch import pandas as pd import numpy as np import pickle

import csv import openpyxl

from ctREFPROP.ctREFPROP import REFPROPFunctionLibrary

Revision 2 : sub-critical flow 추가

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

Reference Table 불러오기

df1 = pd.read_pickle('Pipetable.p')
df2 = pd.read_pickle('Ori_Spring.p')
df3 = pd.read_pickle('Ori_Pilot.p')

Input value load

excel = Dispatch('excel.Application') # win32com.client로 excel call openfiledir = os.getcwd() + '\' # 현재 폴더 위치 5 filename = 'MPB_Ins.xlsx' # 현재 폴더 위치의 해당 파일명의 .xlsx 파일 open

excel.Visible = True

wb = excel.Workbooks.Open(openfiledir + filename) ws = wb.Worksheets("Input")

ws2 = wb.Sheets.Add(Before = None , After = wb.Sheets(wb.Sheets.count))

Composition 설정

Comp_Name = "Nitrogen;Methane;Ethane;Propane;Butane;Isobutane;Pentane;Isopentane"
Comp_Value = [ws.Range("F5").Value, ws.Range("F6").Value, ws.Range("F7").Value, ws.Range("F8").Value, ws.Range("F9").Value, ws.Range("F10").Value, ws.Range("F11").Value, ws.Range("F12").Value]

Input value load

Tag = ws.Range("F2").Value # Tag name, Sheet name 으로 사용.

Fluid = ws.Range("F3").Value

P_set = ws.Range("F16").Value # Set Pressure. P = ws.Range("F17").Value # Relieving Pressure. P1 = ws.Range("F40").Value # Set P + Allowable overpressure + Atmospheric pressure (Inlet Dp 고려??)

F = ws.Range("F25").Value # F factor 는 input 으로 받아야 함. A = ws.Range("F26").Value # Pipe 면적 firefighting = ws.Range("F27").Value # API에 해당하는 인자. firefighting 이 있으면 1 없으면 0 Op_Flow = ws.Range("F37").Value # kg/h operating flowrate (VAPOR 경우, 공급유량 만큼 빼주게 끔)

Tw = ws.Range("F36").Value # K 탄소강판 권장 최대 용기 온도 593 도 Pn = ws.Range("F32").Value # MPa operating pressure Tn = ws.Range("F35").Value # K operating temperature T1 = P/Pn * Tn

Rule = ws.Range("F20").Value # Code 선택 Case = ws.Range("F21").Value # Case 선택 Phase = ws.Range("F22").Value # Phase 선택

Type = ws.Range("F49").Value # Type 선택 Fl_in = ws.Range("F50").Value # Flange rating 선택

P_vent = ws.Range("F42").Value # Vent pressure

단위계 설정

SI = RP.GETENUMdll(0,"SI WITH C").iEnum

Property = "T;CP/CV;Z;M;D;H;"

print (RP.REFPROPdll("AIR","PT","T;H;M;D;T;Z;CP/CV",SI,0,0,0.101325,0,[1,0]).Output[0:Property.count(';')])

print (RP.REFPROPdll("AIR","PT","M;D;H",SI,0,0,0.101325,0,[1,0]).Output[0:Property.count(';')])

Ref. Air 물성치 설정

AIR = np.array(RP.REFPROPdll("AIR","PT",Property,SI,0,0,0.101325,0,[1,0]).Output[0:Property.count(';')])

def Rel_T(a): # if 문으로 liquid, vapor 나누어야 함. Qmass = 0.01

T = RP.REFPROPdll(Comp_Name,"PQmass","XMOLEVAP;XMOLELIQ;HEATVAPZ_P",SI,0,0,a,Qmass,Comp_Value)
Comp_Rel_V = T.Output[0:len(Comp_Value)] # Vapor composition

Comp_Rel_L = T.Output[len(Comp_Value):len(Comp_Value)*2] # Liquid composition

LH = RP.REFPROPdll(Comp_Name,"Pliq","HEATVAPZ_P",SI,0,0,a,0,Comp_Rel_V) # Vapor성분의 lH을 사용으로 변경(R1)
LH = np.array(LH.Output)
Pro_V = RP.REFPROPdll(Comp_Name,"Pvap",Property,SI,0,0,a,0,Comp_Rel_V) # Vapor composition 의 Properties

Pro_L = RP.REFPROPdll(Comp_Name,"Pvap",Property,SI,0,0,a,Q,Comp_Rel_L) # Liquid composition 의 Properties

print (Pro_V.Output[0:Property.count(';')])

Rel_Pro = np.array(Pro_V.Output[0:Property.count(';')])
Rel_Pro = np.append(Rel_Pro,LH[0])
# print(Pro_V)
return (Rel_Pro)

Rel_Pro = Rel_T(P)

def Ins_PRV(a, b): D = math.sqrt(a[1]*(2/(a[1]+1))*((a[1]+1)/(a[1]-1))) D_AIR = math.sqrt(AIR[1](2/(AIR[1]+1))*((AIR[1]+1)/(AIR[1]-1))) W = b (D/D_AIR) math.sqrt((a[3]/AIR[3]) (AIR[0]+273.15)/(a[0]+273.15) * (a[2] / AIR[2]) ) # kg/h return W Ins_nominal = ws.Range("H37").Value W_ins = Ins_PRV(Rel_Pro, Ins_nominal) print(W_ins)

ws_report.Range("F37").Value = W_ins

def Code(a,b,c,d): # Property , b : IGF / API / KOSHA / KGS, c : 'LIQUID' / 'VAPOR', d: Case

if d == 'Fire':    
# 추후 Code 를 dataframe 또는 np.array 로 표로 만들어서 진행하는 방향??
    if b == 'IGF':
        if c == 'LIQUID':
            D = math.sqrt(a[1]*(2/(a[1]+1))**((a[1]+1)/(a[1]-1)))
            D_AIR = math.sqrt(AIR[1]*(2/(AIR[1]+1))**((AIR[1]+1)/(AIR[1]-1)))

            G = 12.4 / (a[6]*D) * math.sqrt (a[2]*(a[0]+273.15)/a[3]) # T는 절대온도로 환산.

            Q = F * G * A**0.82 * 3600 # Nm3/h
            W = Q * AIR[4] * (D/D_AIR) * math.sqrt((a[3]/AIR[3]) * (AIR[0]+273.15)/(a[0]+273.15) * (a[2] / AIR[2]) ) # kg/h
            ## 22. 02. 25 Convert 시 Z 항 추가
            print ("G :" ,G, "Q : ", Q,"W :", W, "Z :", a[2], "Tg", a[0], "Mg", a[3], "D", D)
            return (W)

        elif c == 'VAPOR':
            W = Op_Flow
            return (W)

    elif b == 'API':
        if c == 'LIQUID':        
            if firefighting == 0:
                W = 70900 * F * A**0.82 / a[6] / 1000 * 3600 # kg/h without fire fighting
            elif firefighting == 1:
                W = 43200 * F * A**0.82 / a[6] / 1000 * 3600 # kg/h with fire fighting
            return (W)

        elif c == 'VAPOR':
            W = 0.2772 * math.sqrt(a[3]*P*1000) * (A * (Tw-T1)**1.25 / T1**1.1506)
            return (W)

    elif b == 'KOSHA':
        if c == 'LIQUID':        
            if firefighting == 0:
                W = 4.1868 * 61000 * F * A**0.82 / a[6] # kg/h without fire fighting
            elif firefighting == 1:
                W = 4.1868 * 37000 * F * A**0.82 / a[6] # kg/h with fire fighting        
            return (W)

        elif c == 'VAPOR':
            W = 8.769 * math.sqrt(a[3]*P) * (A * (Tw-T1)**1.25 / T1**1.1506)
            return (W)

    elif b == 'KGS' :
        if c == 'LIQUID':        
            if firefighting == 0:
                W = 4.1868 * 61000 * F * A**0.82 / a[6] # kg/h without fire fighting
            elif firefighting == 1:
                W = 4.1868 * 37140 * F * A**0.82 / a[6] # kg/h with fire fighting          

            return (W)        
        elif c == 'VAPOR':
            W = 0.277 * math.sqrt(a[3]*P*1000) * (A * (Tw-T1)**1.25 / T1**1.1506)
            return (W)
else: # Fire 아닌 경우는 operating condition 으로 배출. 온도만 operating temp.

    W = Op_Flow

S = RP.REFPROPdll(Comp_Name,"PT","S",SI,0,0,Pn,Tn,Comp_Value).Output[0]

Tr = RP.REFPROPdll(Comp_Name,"PS","T",SI,0,0,P,S,Comp_Value).Output[0]

    Rel_Pro[0] = Tn
    return (W)

----------------------- Tank PRV 50% 2개 일때 사용하는 코드 추가 22.04.29

Category = ws.Range("G26").Value # K operating temperature

Rel_Flowrate = np.array(Code(Rel_Pro,Rule,Phase,Case)) ## Required if Category == "Tank": Rel_Flowrate = Rel_Flowrate/2 else: Rel_Flowrate = Rel_Flowrate

-----------------------------

Orifice Sizing

P2 = ws.Range("F42").Value # MPa Kd = ws.Range("F44").Value # C = 0.03948 math.sqrt(Rel_Pro[1] (2/(Rel_Pro[1]+1)) ** ((Rel_Pro[1]+1)/(Rel_Pro[1]-1)))

Kb : initial 1 (for conventional, pilot type), calculated value 2()

Kb = np.array([1, (1/(17.9 C) math.sqrt((Rel_Pro[1]/(Rel_Pro[1]-1)) * (P2/P)(2/Rel_Pro[1]) - (P2/P)((Rel_Pro[1]+1)/Rel_Pro[1]) ))])

kc 는 1을 보통 쓰나, critical flow 가 나뉘면서 metrix 화되어, 공통 인자는 kc를 list 화 시킴.

------------------------kc 가 왜 2개로 되어있는지 모르겠음. 하나로 수정해서 Req_ori 도 list 가 아닌 float 으로만 남김. (22.04.29)

Kc = ws.Range("F45").Value

Kc = np.array([ws.Range("F45").Value, 1]) # Combination correction factor of installations with a rupture disk, No rupture disk = 1

Critical flow check

P_criflow = P1000 (2/(Rel_Pro[1]+1))*(Rel_Pro[1]/(Rel_Pro[1]-1)) if P_criflow > P21000 : # Critical flow 에서는 P1 이 relieving pressure. Inlet pressure 가야하나? Check_cri = "Critical flow"

Required Orifice size (kb 에 따라 2개 list 로 정의함.)

Req_Ori = Rel_Flowrate / (C * Kd * P*1000 * Kb[0] * Kc) * math.sqrt((Rel_Pro[0]+273.15)*Rel_Pro[2]/Rel_Pro[3]) # critical flow

else: # Sub critical flow 에서는 P1 을 Valve Inlet pressure 로 정의. Inlet pressure drop 도 고려해서 orifice size 크게 가져간다. Check_cri = "Sub-Critical flow" r = P2/P1 F2 = math.sqrt( (Rel_Pro[1]/(Rel_Pro[1]-1))* r*(2/Rel_Pro[1]) (1- r*((Rel_Pro[1]-1)/Rel_Pro[1]))/(1-r) ) Req_Ori = 17.9 Rel_Flowrate / (F2 Kd Kc) math.sqrt((Rel_Pro[0]+273.15)Rel_Pro[2]/ (Rel_Pro[3] P11000(P11000 - P21000))) # Sub-critical flow
Req_Ori_ref = Rel_Flowrate / (C Kd P1000 Kb[1] Kc) math.sqrt((Rel_Pro[0]+273.15)Rel_Pro[2]/Rel_Pro[3]) print("r :", r, "F2 :", F2,"P1 :", P1,"P2 :", P2)

Required Orifice size (kb 에 따라 2개 list 로 정의함.)

Req_Ori = Rel_Flowrate / (C Kd P1000 Kb Kc) math.sqrt((Rel_Pro[0]+273.15)*Rel_Pro[2]/Rel_Pro[3]) # critical flow

r = P2/P

F2 = math.sqrt( (Rel_Pro[1]/(Rel_Pro[2]-1))* r*(2/Rel_Pro[2]) (1- r**((Rel_Pro[2]-1)/Rel_Pro[2]))/(1-r) )

Req_Ori = 17.9 Rel_Flowrate / (F2 Kd Kc) math.sqrt((Rel_Pro[0]+273.15)Rel_Pro[2]/ (Rel_Pro[3] P1000(P1000 - P21000))) # Sub-critical flow

API 520의 Orifice Size, Designation 구하는 함수

def Ori(a): # a : Req_Ori[0] Size = np.array([['D', 0.000001, 70.97], ['E', 70.97, 126.45], ['F', 126.45, 198.06], ['G', 198.06, 324.52], ['H', 324.52, 506.45], ['J', 506.45, 830.32], ['K', 830.32, 1185.8], ['L', 1185.8, 1840.64], ['M', 1840.64, 2322.58], ['N', 2322.58, 2799.99], ['P', 2799.99, 4116.12], ['Q', 4116.12, 7129.02], ['R', 7129.02, 10322.56], ['T', 10322.56, 16774.16]] )

Ori_A = np.array

for i in Size:
    # print (a,Ori_A)
    if a < float(Size[0,1]): # 최소값 error항
        Ori_A = "Out of API range(under)"
    elif a > float(Size[13,2]): # 최대값 error항
        Ori_A = "Out of API range(over)"
    elif a > float(i[1]): # Orifice size 찾는 항
        Ori_A = (i[0], float(i[2]))
    # print(i, Ori_A)

return Ori_A # Ori_A 에는 designation, Orifice area 정보 포함.

Required Orifice size (kb 에 따라 2개 list 로 정의함.)

대부분 0을 쓰나, subcritical 경우에는 vendor 에 문의, API520표에 보면 bellows 에서 10% overpressure 일대는 표 참조.

------------------- Cri, Sub-cri 로 이미 계산되어 kb_f 제외 22.04.29

kb_f = int(ws.Range("F46").Value)

Ori_A = Ori(Req_Ori)

def Flange(a,b): # a : Spring / Pilot, b : Flange size

print (a,b)

if a == "Bellows" or a == "Conventional":
    Ori_Table = df2.loc[Ori(Req_Ori)[0]] # Orifice designation 에 맞는 table 찾기
    df_Orifice = Ori_Table.loc[Ori_Table["Flange_Inlet"] == b, :]

print (Orifice)

    return df_Orifice

elif a == "Pilot Operated":
    Ori_Table = df3.loc[Ori(Req_Ori)[0]]
    # Ori_Table = df3.loc[Ori(Req_Ori)[0]]
    df_Orifice = Ori_Table.loc[Ori_Table["Flange_Inlet"] == b, :]   
    print(Ori_Table, df_Orifice)

print (Orifice)

    return df_Orifice

kc는 critical flow list 맞추기 위해 사용. 0,1 모두 1의 값을 가짐

if Ori_A == "Out of API range(under)" or Ori_A =="Out of API range(over)": data_flange = [["-","-","-","-"]] df_Orifice = pd.DataFrame(data_flange, columns =['Flange_Inlet','Orifice_Inlet','OrificeOutlet','Flange Outlet']) Rated_flowrate = "-" df_Orifice["Required Capacity (kg/h)"] = Rel_Flowrate df_Orifice["Orifice Area (m2)"] = Ori_A df_Orifice["Code"] = Rule df_Orifice["Case"] = Case df_Orifice["Phase"] = Phase

else:

---------------------------------------- Cri, Sub-cri 에 따라 kb 값 사용

if Check_cri =='Critical flow':
    kb_f = 0
else:
    kb_f = 1

------------------------------------------

Rated_flowrate = Ori_A[1] * C * Kd * P*1000 * Kb[kb_f] * Kc * math.sqrt(Rel_Pro[3]/(Rel_Pro[0]+273.15)/Rel_Pro[2]) 
                          #(C * Kd * P*1000 * Kb[0] * Kc) * math.sqrt((Rel_Pro[0]+273.15)*Rel_Pro[2]/Rel_Pro[3])
df_Orifice = Flange(Type,Fl_in)
df_Orifice["Required Capacity (kg/h)"] = Rel_Flowrate
df_Orifice["Orifice Area (m2)"] = Ori_A[1]
df_Orifice["Code"] = Rule
df_Orifice["Case"] = Case
df_Orifice["Phase"] = Phase

df_Orifice = Flange(Type,Fl_in)

df_Orifice["Required Capacity (kg/h)"] = Rel_Flowrate

df_Orifice["Orifice Area (m2)"] = Ori_A[1]

df_Orifice["Code"] = Rule

df_Orifice["Case"] = Case

df_Orifice["Phase"] = Phase

print (df_Orifice)

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' Report 작성 '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

wb.Sheets("Format").Copy(Before = None, After = wb.Sheets("Format")) wb.Activesheet.Name = Tag ws_report = wb.Worksheets(Tag)

ws_report.Range("L6").Value = Tag ws_report.Range("L7").Value = Type ws_report.Range("L8").Value = Rule ws_report.Range("AA6").Value = Case ws_report.Range("AA8").Value = Fl_in if Ori_A == "Out of API range(under)" or Ori_A =="Out of API range(over)": ws_report.Range("AD8").Value = dfOrifice['Flange Outlet'] ws_report.Range("L13").Value = df_Orifice['Orifice Area (m2)'] ws_report.Range("AA9").Value = df_Orifice['Required Capacity (kg/h)'] ws_report.Range("AA12").Value = Ori_A else: ws_report.Range("AD8").Value = float(dfOrifice['Flange Outlet']) ws_report.Range("L13").Value = float(df_Orifice['Orifice Area (m2)']) ws_report.Range("AA9").Value = float(df_Orifice['Required Capacity (kg/h)']) ws_report.Range("AA12").Value = Ori_A[0] ws_report.Range("L9").Value = F ws_report.Range("L10").Value = Rel_Pro[6] ws_report.Range("L11").Value = A ws_report.Range("L12").Value = Req_Ori ws_report.Range("AA10").Value = Rated_flowrate ws_report.Range("AA11").Value = Rel_Pro[0]

ws_report.Range("AA13").Value = df_Orifice['Orifice_Inlet'] ws_report.Range("AA14").Value = P_vent1000 ws_report.Range("AD13").Value = df_Orifice['Orifice_Outlet'] ws_report.Range("L14").Value = Fluid ws_report.Range("L15").Value = Rel_Pro[3] ws_report.Range("L16").Value = Rel_Pro[1] ws_report.Range("L17").Value = Rel_Pro[2] ws_report.Range("L19").Value = P_set ws_report.Range("L21").Value = P1000

ws_report.Range("AA16").Value = Kb[0]

ws_report.Range("AA20").Value = Kc ws_report.Range("AA21").Value = Kd ws_report.Range("AA22").Value = C ws_report.Range("AA16").Value = Check_cri

if Check_cri =="Sub-Critical flow": ws_report.Range("AB34").Value = P11000 - 101.325 else: ws_report.Range("AB34").Value = P1000 - 101.325

Composition write

C = RP.REFPROPdll(Comp_Name,"PQmass","XMOLEVAP;XMOLELIQ;HEATVAPZ_P",SI,0,0,P,0.01,Comp_Value) Composition = C.Output[0:len(Comp_Value)] # Vapor composition

ws_report.Range("E44").Value = Composition[0] ws_report.Range("E45").Value = Composition[1] ws_report.Range("E46").Value = Composition[2] ws_report.Range("E47").Value = Composition[3] ws_report.Range("E48").Value = Composition[4] ws_report.Range("E49").Value = Composition[5] ws_report.Range("E50").Value = Composition[6] ws_report.Range("E51").Value = Composition[7]

df1 = pd.read_pickle('Pipetable.p')

def pipe(a,b): # ND, 'Sch'

print ("ND".center(8),"Schedule".center(8),"Inner_Dia.".center(12), "Out_Dia.".center(12), "Thickness".center(12))

print (a, b, df1.loc[a,"Out_Dia"]-2*df1.loc[a,b], df1.loc[a,"Out_Dia"], df1.loc[a,b])

pipe(50,'10s')

#

JungHulk commented 2 years ago

Need to add ammonia and carbon dioxide mode.