気液二相流における配管圧力損失計算ソフト

import numpy as np
import pandas as pd 
import glob
import math

data_sheet = pd.read_excel("steam data input.xlsx",index_col=0)　

#データシートを読み込む
data_sheet = pd.read_excel("steam data input.xlsx",index_col=0)
ds = data_sheet
ds.columns = ["result","unit"]

piping_data = pd.read_excel("steam data input.xlsx",sheet_name='配管内径　厚さ',header=2,index_col=0)
pp =  piping_data
pp.columns = ["外径","sh5_thk","sh5_di","sh10_thk","sh10_di","sh20_thk","sh20_di","sh40_thk","sh40_di","0","sh80_thk","sh80_di","sh120_thk","sh120_di","sh160_thk","sh160_di"]

P = ds.at['圧力','result']
flowrate_steam = ds.at['蒸気流量','result']
flowrate_hotwater = ds.at['熱水流量','result']
length_piping = ds.at['配管長','result']
L = length_piping
piping = ds.at['配管径','result']

di = pp.at[str(piping),'sh10_di']
di = di*0.001

#流量
Q = (flowrate_steam + flowrate_hotwater)*1000
Q = Q/3600

#乾き度
X = flowrate_steam/(flowrate_steam + flowrate_hotwater)


"""
蒸気表をあるwebサイトからseleniumuを用いて、引用しているためここでは省略します

"""

Temp = steam_data[0].text
latent_heat = steam_data[1].text
enthalpy_steam = steam_data[2].text
enthalpy_hotwater = steam_data[3].text
specific_volume_steam = steam_data[4].text
specific_volume_hotwater = steam_data[5].text

temp = float(Temp)
enthalpy_steam = float(enthalpy_steam)
enthalpy_hotwater = float(enthalpy_hotwater)
latent_heat = float(latent_heat)

#密度　kg/m3
density_G = 1/float(specific_volume_steam)
density_L = 1/float(specific_volume_hotwater)

#蒸気粘性関数
"""
viscosity_steam = pd.read_excel("steam data input.xlsx",sheet_name='蒸気粘性',header=1,index_col=0)

viscosity_steam = viscosity_steam.query('0 < MPaG < 1.5')
x = viscosity_steam['MPaG']
y = viscosity_steam['mPa・s']

cf1 = np.polyfit(x,y,1)
cf2 = np.polyfit(x,y,2)

func1 = np.poly1d(cf1)
print(func1)
fig, ay = plt.subplots()
ay.invert_yaxis()
plt.plot(x, func1(x),label='d=1')

func2 = np.poly1d(cf2)
print(func2)
plt.plot(x, func2(x),label='d=2')

plt.scatter(x,y)
plt.grid()
"""
myu_G = -0.00127*P**2 + 0.003889*P + 0.01263

#熱水粘性関数
""""
viscosity_hotwater = pd.read_excel("steam data input.xlsx",sheet_name='熱水粘性',header=1)
viscosity_hotwater1 = viscosity_hotwater.query('0 <Temp< 100')
x = viscosity_hotwater1['Temp']
y = viscosity_hotwater1['mPa・s']

cf1 = np.polyfit(x,y,1)
cf2 = np.polyfit(x,y,2)


func1 = np.poly1d(cf1)
print(func1)
fig, ay = plt.subplots()
ay.invert_yaxis()
plt.plot(x, func1(x),label='d=1')

func2 = np.poly1d(cf2)
print(func2)
plt.plot(x, func2(x),label='d=2')

plt.scatter(x,y)
plt.grid()
"""

t= temp
if temp > 100:
    myu_L = 4.665e-06*t**2 - 0.002739*t + 0.4946
else:
    myu_L = 0.0001644*t**2 - 0.02834*t + 1.546

#二相流圧力損失計算　
#Homogeneous Flow Model Applied to Intube Flow
#The homogeneous void fraction
E = 1/(1+((1-X)/X)*(density_G/density_L))
print(E)

#均質密度
RoH = density_L*(1-E)+density_G*E
print(RoH)

#均質粘性
myu_tp = X*myu_G + (1-X)*myu_L
print(myu_tp)

#レイノルズ数
S = math.pi*(di/2)**2
m = Q/S
Re =m*di/myu_tp #1mの場合、Ｌが長くなると乱流になりやすい。Re増
print('S = ' + str(S)+' m2')
print('m = ' + str(m)+ ' kg/m s')
print('Re = '+ str(Re))

if Re > 2300:
    statusH = '乱流'
    print(statusH)
    f_tp = 0.079/(Re**0.25)
    #f_tp = 0.046/(Re**0.20)
elif Re<2300:
    statusH = '層流'
    print(statusH)
    f_tp = 16/Re

#f_tp1 = 0.079/(Re**0.25)
#f_tp2 = 0.046/(Re**0.20)    

print('f_tp = ' + str(f_tp)) #摩擦係数

#圧力損失
delta_P = 2*f_tp*(m**2)/(di*RoH)
print(str(delta_P)+ ' N/m2')

delta_P =delta_P/1000
print(str(delta_P)+ ' kPa')

delta_P =delta_P/1000
print(str(delta_P)+ ' MPa')

#見かけレイノルズ数

QG = flowrate_steam*1000/3600
QL = flowrate_hotwater*1000/3600

mG = QG/S
mL = QL/S


ReG =mG*di/myu_G
ReL =mL*di/myu_L

print('Re(気体) = '+ str(ReG))
print('Re(液体) = '+ str(ReL))

if ReG > 1500:
    if ReL >1500:
        status = '気相、液相ともに乱流、t-t'
        print(status)
    else:
        status = '気相－乱流、液相－層流、t-v'
        print(status)
else:
    if ReL >1500:
        status = '気相－層流、液相－乱流、v-t'
        print(status)
    else:
        status = '気相－層流、液相－層流、v-v'
        print(status)

#摩擦損失係数
if ReG <= 1000:
    f_G = 64/ReG
else:
    f_G = 0.079*(ReG**-0.25)

if ReL <= 1000:
    f_L = 64/ReL
else:
    f_L = 0.079*(ReL**-0.25)
    
print(f_G)
print(f_L)

#見かけの圧力損失
dP_G = 2*f_G*(mG**2)/(di*density_G)
dP_L = 2*f_L*(mL**2)/(di*density_L)

XX = (dP_L/dP_G)**0.5

print(dP_G)
print(dP_L)
print(XX)

#グラフ読み取り　アナログ
fai_G=1.8
fai_L=3.1

dP1 = (fai_G**2)*dP_G
dP1 = dP1*(10**-6)

print('dP1 = ' +str(dP1) + ' MPaG/m')

dP2 = (fai_L**2)*dP_L
dP2 = dP2*(10**-6)
print('dP2 = ' +str(dP2) + ' MPaG/m')

dP = (dP1+dP2)/2
print('dP = ' +str(dP) + ' MPaG/m')

import openpyxl as px
import pprint

wb = px.load_workbook('two phrase pressure drop output.xlsx')
ws = wb.active
sheet = wb['two phrase pressure drop']



#ある列に順に書き込む関数
def write_list_1d(sheet, l_1d, start_row, start_col):
    for i, row in enumerate(l_1d):
        sheet.cell(row=start_row + i,column=start_col, value=l_1d[i])
        
        
l_1d = [P,temp, flowrate_steam, flowrate_hotwater, L , piping]
write_list_1d(sheet, l_1d, 2, 3)

#蒸気表の記入
l_1d1 = [temp, 
         latent_heat, enthalpy_steam, enthalpy_hotwater, 
         density_G,density_L,
         X,myu_G,myu_L]
write_list_1d(sheet, l_1d1, 2, 9)

#均質モデル　パラメータ
l_1d2 = [E,RoH, myu_tp, m]
write_list_1d(sheet, l_1d2, 13, 9)

#均質モデル　結果
l_1d3 = [Re,statusH, f_tp, delta_P]
write_list_1d(sheet, l_1d3, 13, 14)

#分離モデル　パラメータ
l_1d4 = [mG,mL]
write_list_1d(sheet, l_1d4, 19, 9)

#分離モデル　結果
l_1d5 = [ReG, ReL, status, f_G,f_L,XX, fai_G, fai_L ,dP] 
write_list_1d(sheet, l_1d5, 19, 14)
        
wb.save('two phrase pressure drop output.xlsx')