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

SJ_PRV #4

Open
JungHulk opened this issue Sep 7, 2022 · 1 comment
Open

SJ_PRV #4

JungHulk opened this issue Sep 7, 2022 · 1 comment

Comments

@JungHulk
Copy link
Owner

JungHulk commented Sep 7, 2022

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 * P1000 * 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 * P*1000 * 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 * P*1000 * 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','Orifice_Outlet','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 * P1000 * Kb[kb_f] * Kc * math.sqrt(Rel_Pro[3]/(Rel_Pro[0]+273.15)/Rel_Pro[2])
#(C * Kd * P1000 * 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 = df_Orifice['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(df_Orifice['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
Copy link
Owner Author

JungHulk commented Sep 7, 2022

Need to add ammonia and carbon dioxide mode.

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