File:
```python
# -*- coding: utf-8 -*-
"""
Created on Sun, April 17 2022
This script is written to analytically and semi-analytical-emperically fitting to solve overall kGAv for spray system that absorbs CO2 from the surrounding air/flue gas
@author: Awan
"""
import math
import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import *
import pandas as pd
from Input_Params import *
from Side_Functions import *
plt.close('all')
WritetoFolder='Results'
# def KGAv_Theo(Q_l_flux, Q_g_flux, P_CO2_in, T, CO2_loading):
# hm, U, d = hm_calc(Q_l_flux, Q_g_flux, P_CO2_in, T, CO2_loading)
# A_1d = 4*3.14*d*d/4
# Vol_1d = (4/3)*3.14*(d**3)/8
# Q_L = Q_l_flux*A_D/3600 # Liquid flow rate [m3/s] when Q_l_flux is in m3/m2-h
# Q_G =Q_g_flux*A_D/3600 # Gas Flow Rate [m3/s] (For aroonwilas paper, the flow rate is assumed total gas flow rate for now as info is not clerarly mentioned in paper)
# Q_G_mole = 3600*(Q_G*1000/(22.4*1000))*((T/298)/(P_tot/101.325))/(A_D) #Gas Flow Rate in kmol/m2-hr
# X2_theo = X1*math.exp(-hm*Q_L*A_1d*H/(Q_G*Vol_1d*U))
# Y2_theo = X2_theo/(1-X2_theo)
# KGAv_theo = (Q_G_mole*(1-X1)/(P_tot*H))*(math.log(Y1/Y2_theo) + Y1 - Y2_theo)
# cap_eff = (X1-X2_theo)*100/X1 ### Change the definition as per paper that it is being validated to, if needed
# return KGAv_theo, cap_eff
X = np.zeros(2) #X[0] is loading and X[1] is X_CO2
# def loading_func(z, CO2_loading):
# hm, U, d = hm_calc(Q_l_flux, Q_g_flux, P_CO2_in, T, CO2_loading)
# AB = (P_tot/(C_init*Rbar*T))*(6/(d*U))*hm
# return AB*P_CO2_in/P_tot
# def X_CO2_func(z,X_CO2):
# hm, U, d = hm_calc(Q_l_flux, Q_g_flux, P_CO2_in, T, CO2_loading)
# B = (Q_l_flux/Q_g_flux)*(6/(d*U))*hm
# return B*X_CO2
def X_func(z, X):
hm, d = hm_calc(Q_l_flux, Q_g_flux, X[1]*P_tot, T, X[0]/(1-X[0]), U)
AB = (P_tot/(C_init*Rbar*T))*(6/(d*U))*hm
B = (Q_l_flux/Q_g_flux)*(6/(d*U))*hm
print(AB*X[1], B*X[1])
return AB*X[1], B*X[1]
npoints = len(x_exp)
Inputs = np.zeros([len(x_exp),4])
KGAv_Theo_all = np.zeros(npoints)
cap_eff_all = np.zeros(npoints)
KGAv_Fit_all = np.zeros(npoints)
cap_eff_Fit_all = np.zeros(npoints)
# def X2_calc_integrand():#, Q_l_flux, Q_g_flux, P_CO2_in, T):
# hm, U, d = hm_calc(Q_l_flux, Q_g_flux, P_CO2_in, T, CO2_loading_bottomUp)
# A_1d = 4*3.14*d*d/4
# Vol_1d = (4/3)*3.14*(d**3)/8
# Q_L = Q_l_flux*A_D/3600 # Liquid flow rate [m3/s] when Q_l_flux is in m3/m2-h
# Q_G =Q_g_flux*A_D/3600 # Gas Flow Rate [m3/s] (For aroonwilas paper, the flow rate is assumed total gas flow rate for now as info is not clerarly mentioned in paper)
# return -hm*Q_L*A_1d/(Q_G*Vol_1d*U)
for i in range(npoints):
if Var == 'Liquid Flow Rate':
Q_l_flux = x_exp[i]
Inputs[i] = [Q_l_flux, Q_g_flux, P_CO2_in, T]
if Var == 'Gas Flow Rate':
Q_g_flux = x_exp[i]
Inputs[i] = [Q_l_flux, Q_g_flux, P_CO2_in, T]
U = U_calc2(Q_l_flux, Q_g_flux, T)
print('Starting X_sol')
# loading_sol = scipy.integrate.solve_ivp(loading_func, [0, H], [0])
# X_CO2_sol = scipy.integrate.solve_ivp(X_CO2_func, [0, H], [0.01])
X_sol = scipy.integrate.solve_ivp(X_func, (0, H), [0, 0.1])
# CO2_loading_bottomUp = loading_sol.y[::-1]
# z_array_bottomUp = loading_sol.t[::-1]
# z_array_bottomUp[:] = H-z_array_bottomUp[:]
# hm, U, d = hm_calc(Q_l_flux, Q_g_flux, P_CO2_in, T, CO2_loading_bottomUp)
# ln_X2_array = X2_calc_integrand()
# Integrand = scipy.integrate.simps(ln_X2_array, z_array_bottomUp)
# print (Integrand)
# X2 = X1*math.exp(Integrand)
#lnX2_by_X1 = scipy.integrate.quad(X2_calc_integrand, 0, H)
#print(X2)
# for j in range(len(loading_sol.y)):
# hm_array =
# KGAv_theo, cap_eff = KGAv_Theo(Q_l_flux, Q_g_flux, P_CO2_in, T, CO2_loading)
# KGAv_Theo_all[i] = KGAv_theo
# cap_eff_all[i] = cap_eff
# def KGAv_Fit(Inputs, K, c):
# for i in range(npoints):
# [Q_l_flux, Q_g_flux, P_CO2_in, T] = Inputs[i]
# hm, U, d = hm_calc(Q_l_flux, Q_g_flux, P_CO2_in, T)
# A_1d = 4*3.14*d*d/4
# Vol_1d = (4/3)*3.14*(d**3)/8
# Q_L = Q_l_flux*A_D/3600 # Liquid flow rate [m3/s] when Q_l_flux is in m3/m2-h
# Q_G =Q_g_flux*A_D/3600 # Gas Flow Rate [m3/s] (For aroonwilas paper, the flow rate is assumed total gas flow rate for now as info is not clerarly mentioned in paper)
# Q_G_mole = 3600*(Q_G*1000/(22.4*1000))*((T/298)/(P_tot/101.325))/(A_D) #Gas Flow Rate in kmol/m3-hr
# X2_Fit = X1*math.exp(-hm*Q_L*A_1d*H*K/(Q_G*Vol_1d*U))
# Y2_Fit = X2_Fit/(1-X2_Fit)
# KGAv_Fit = (Q_G_mole*(1-X1)/(P_tot*H))*(math.log(Y1/Y2_Fit) + Y1 - Y2_Fit) + c
# KGAv_Fit_all[i] = KGAv_Fit
# cap_eff_Fit = (X1-X2_Fit)*100/X1 ### Change the definition as per paper that it is being validated to, if needed
# cap_eff_Fit_all[i] = cap_eff_Fit
# return KGAv_Fit_all
# params, extras = curve_fit(KGAv_Fit, Inputs, y_exp)
# K = params[0]
# c = params[1]
# print(K,c)
# KGAv_Fit_all = KGAv_Fit(Inputs, K, c)
plt.figure()
plt.plot(x_exp, KGAv_Theo_all, 'go-', label = 'Analytical Model')
# plt.plot(x_exp, KGAv_Fit_all, 'ro-', label = 'Model with Fitting')
plt.plot(x_exp, y_exp, 'bo-', label = 'Expt Data')
plt.legend()
# if Var == 'Liquid Flow Rate':
# Q_l_flux_array = np.zeros(len(x_exp))
# Q_l_flux_array[:] = x_exp[:]
# Q_g_flux = 382 #m3/m2-h
# Q_G =Q_g_flux*A_D/3600 #Gas Flow Rate [m3/s] (Assumed total gas flow rate for now as info is not clerarly mentioned in paper)
# Q_G_mole = 3600*(Q_G*1000/(22.4*1000))*(((273+T)/298)/(P_tot/101.325))/(A_D) #Gas Flow Rate in kmol/m3-hr
# hm_array, U_array, d_array = hm_calc_QLarray(Q_l_flux_array, Q_g_flux, A_D, K_factor, d_orifice, P_CO2_in, P_tot, T, nozzle_type)
# def KGAv_Theo_for_Ql(Q_l_flux_array, hm_array, U_array, d_array):
# Q_L_array = np.zeros(len(Q_l_flux_array))
# X2_theo = np.zeros(len(Q_l_flux_array))
# Y2_theo = np.zeros(len(Q_l_flux_array))
# KGAv_theo = np.zeros(len(Q_l_flux_array))
# for i in range(len(Q_l_flux_array)):
# Q_L_array[i] = Q_l_flux_array[i]*A_D/3600
# d = d_array[i]
# A_1d = 4*3.14*d*d/4
# Vol_1d = (4/3)*3.14*(d**3)/8
# X2_theo[i] = X1*math.exp(-hm_array[i]*Q_L_array[i]*A_1d*H/(Q_G*Vol_1d*U_array[i]))
# #X2_calc[i] = X1*math.exp(-c*H*((Q_L_array[i]*d_orifice)**.5)/Q_G)
# Y2_theo[i] = X2_theo[i]/(1-X2_theo[i])
# KGAv_theo[i] = (Q_G_mole*(1-X1)/(P_tot*H))*(math.log(Y1/Y2_theo[i])+Y1-Y2_theo[i])
# return KGAv_theo
# if Var == 'Gas Flow Rate':
# Q_g_flux_array = np.zeros(len(x2_exp))
# Q_g_flux_array[:] = x2_exp[:]
# Q_L = Q_l_flux*A_D/3600
# hm_array, U_array, d_array = hm_calc_QGarray(Q_g_flux_array, Q_l_flux, A_D, K_factor, d_orifice, P_CO2_in, P_tot, T, nozzle_type)
# def KGAv_Theo_for_Qg(Q_g_flux_array, hm_array, U_array, d_array):
# Q_G_array = np.zeros(len(Q_g_flux_array))
# Q_G_mole_array = np.zeros(len(Q_g_flux_array))
# X2_theo = np.zeros(len(Q_g_flux_array))
# Y2_theo = np.zeros(len(Q_g_flux_array))
# KGAv_theo = np.zeros(len(Q_g_flux_array))
# for i in range(len(Q_g_flux_array)):
# Q_G_array[i] = Q_g_flux_array[i]*A_D/3600
# Q_G_mole_array[i] = 3600*(Q_G_array[i]*1000/(22.4*1000))*(((273+T)/298)/(P_tot/101.325))/(A_D)
# d = d_array[i]
# A_1d = 4*3.14*d*d/4
# Vol_1d = (4/3)*3.14*(d**3)/8
# X2_theo[i] = X1*math.exp(-hm_array[i]*Q_L*A_1d*H/(Q_G_array[i]*Vol_1d*U_array[i]))
# #X2_calc[i] = X1*math.exp(-c*H*((Q_L_array[i]*d_orifice)**.5)/Q_G)
# Y2_theo[i] = X2_theo[i]/(1-X2_theo[i])
# KGAv_theo[i] = (Q_G_mole_array[i]*(1-X1)/(P_tot*H))*(math.log(Y1/Y2_theo[i])+Y1-Y2_theo[i])
# return KGAv_theo
# def KGAv_Calculator_for_Ql(Q_l_flux_array, c, K):
# Q_L_array = np.zeros(len(Q_l_flux_array))
# X2_calc = np.zeros(len(Q_l_flux_array))
# Y2_calc = np.zeros(len(Q_l_flux_array))
# KGAv_calc = np.zeros(len(Q_l_flux_array))
# for i in range(len(Q_l_flux_array)):
# Q_L_array[i] = Q_l_flux_array[i]*A_D/3600
# X2_calc[i] = X1*math.exp(-c*hm_array[i]*H*((Q_L_array[i]*d_orifice)**.5)/Q_G)
# #X2_calc[i] = X1*math.exp(-c*H*((Q_L_array[i]*d_orifice)**.5)/Q_G)
# Y2_calc[i] = X2_calc[i]/(1-X2_calc[i])
# #KGAv_calc[i] = (Q_G_mole*(1-X1)/(P_tot*H))*(math.log(Y1/Y2_calc[i])+Y1-Y2_calc[i])
# KGAv_calc[i] = (Q_G_mole*(1-X1)/(P_tot*H))*(math.log(Y1/Y2_calc[i])+Y1-Y2_calc[i]) + K
# return KGAv_calc
# #def KGAv_Calculator_for_Qg(Q_g_flux_array, c, K):
# def KGAv_Calculator_for_Qg(Q_g_flux_array, c, K, c2):
# Q_G_array = np.zeros(len(Q_g_flux_array))
# X2_calc = np.zeros(len(Q_g_flux_array))
# Y2_calc = np.zeros(len(Q_g_flux_array))
# KGAv_calc = np.zeros(len(Q_g_flux_array))
# for i in range(len(Q_g_flux_array)):
# Q_G_array[i] = Q_g_flux_array[i]*A_D/3600
# X2_calc[i] = X1*math.exp(-c*hm_array[i]*H*((Q_L*d_orifice)**.5)/Q_G_array[i]**c2)
# #X2_calc[i] = X1*math.exp(-c*6*hm_array[i]*H*Q_L/(d_array[i]*Q_G_array[i]*U_array[i]))
# Y2_calc[i] = X2_calc[i]/(1-X2_calc[i])
# print(Y2_calc[i])
# Q_G_mole = 3600*(Q_G_array[i]*1000/(22.4*1000))*(((273+T)/298)/(P_tot/101.325))/(A_D)
# #KGAv_calc[i] = (Q_G_mole*(1-X1)/(P_tot*H))*(math.log(Y1/Y2_calc[i])+Y1-Y2_calc[i])
# KGAv_calc[i] = (Q_G_mole*(1-X1)/(P_tot*H))*(math.log(Y1/Y2_calc[i])+Y1-Y2_calc[i]) + K
# return KGAv_calc
# if Var == 'Liquid Flow Rate':
# params, extras = curve_fit(KGAv_Calculator_for_Ql, x_exp, y_exp)
# c = params[0]
# K = params[1]
# plt.plot(x_exp, KGAv_Theo_for_Ql(Q_l_flux_array, hm_array, U_array, d_array), 'go-', label = 'Analytical Model')
# plt.plot(x_exp, KGAv_Calculator_for_Ql(Q_l_flux_array, c, K), 'ro-', label = 'Semi-Analytical-Empirical Fit')
# plt.plot(x_exp, y_exp, 'bo-', label = 'Experimental data')
# plt.xlabel('Liquid flow rate (m3/m2-h)')
# if Var == 'Gas Flow Rate':
# params, extras = curve_fit(KGAv_Calculator_for_Qg, x2_exp, y_exp)
# c = params[0]
# K = params[1]
# c2 = params[2]
# #plt.plot(x2_exp, KGAv_Theo_for_Qg(Q_g_flux_array, hm_array, U_array, d_array), 'g', label = 'Theoretical Model')
# plt.plot(x2_exp, KGAv_Calculator_for_Qg(Q_g_flux_array, c, K, c2), 'ro-', label = 'Analytical Fit')
# #plt.plot(x2_exp, KGAv_Calculator_for_Qg(Q_g_flux_array, c, K), 'ro', label = 'Analytical Fit')
# plt.plot(x2_exp, y_exp, 'bo-', label = 'Experimental data')
# plt.xlabel('Gas flow rate (m3/m2-h)')
# plt.ylabel('KgAe (kmol/m3-h-kPa)')
# plt.legend(loc='best')
# #plt.title('Kuntz fit 2f and 2c')
# #plt.savefig()#'%s/Concentration_timehistory.pdf'%WritetoFolder, bbox_inches='tight')
# if nozzle_type == 'P-40':
# c_P40 = c
# K_P40 = K
# if nozzle_type == 'P-28':
# c_P28 = c
# K_P28 = K
# print(c, K)
######plt.plot(x_exp, KGAv_Calculator_for_Ql(x_exp, c), 'r')
# plt.plot(x_exp, KGAv_Calculator_for_Ql(x_exp, c, K), 'ro', label = 'Analytical Fit')
# #####plt.plot(x_exp, KGAv_Calculator_for_Ql(x_exp, 504.1583049223172, -1.6742542161512164), 'g')
# #####plt.plot(x_exp, KGAv_Calculator_for_Ql(x_exp, 246.1363698486537, -0.17978853997655872), 'g')
# ######plt.plot(x_exp, KGAv_Calculator_for_Ql(x_exp, 1737430449.8144524, 1.39), 'g')
# plt.plot(x_exp, y_exp, 'bo', label = 'Experimental data')
#plt.plot(x2_exp, KGAv_Calculator_for_Qg(x2_exp, c,c2, K), 'ro', label = 'Analytical Fit')
#plt.plot(x2_exp, y_exp, 'bo', label = 'Experimental data')
#
# r_probe=R/8
# t_probe=1
# C_probe=DropletDiffusion_Functions.calc_probe(r_probe,t_probe,R,Rbar,kg,Bi,Henry,DCO2,Cinit,PCO2)
# tstart=1
# tfinal=100
# Ntime=20
# tsample=np.linspace(tstart,tfinal,Ntime)
# C_all_Rover4=np.zeros((Ntime,))
# C_all_Rover2=np.zeros((Ntime,))
# C_all_Rover1=np.zeros((Ntime,))
# for i in range(0,Ntime):
# C_all_Rover4[i]=DropletDiffusion_Functions.calc_probe(R/4,tsample[i],R,Rbar,kg,Bi,Henry,DCO2,Cinit,PCO2)
# C_all_Rover2[i]=DropletDiffusion_Functions.calc_probe(R/2,tsample[i],R,Rbar,kg,Bi,Henry,DCO2,Cinit,PCO2)
# C_all_Rover1[i]=DropletDiffusion_Functions.calc_probe(R/1,tsample[i],R,Rbar,kg,Bi,Henry,DCO2,Cinit,PCO2)
# Ngrid=20
# rstart=R/Ngrid
# rsample=np.linspace(rstart,R,Ngrid)
# C_all_t10=np.zeros((Ngrid,))
# C_all_t20=np.zeros((Ngrid,))
# C_all_t30=np.zeros((Ngrid,))
# for i in range(0,Ngrid):
# C_all_t10[i]=DropletDiffusion_Functions.calc_probe(rsample[i],10,R,Rbar,kg,Bi,Henry,DCO2,Cinit,PCO2)
# C_all_t20[i]=DropletDiffusion_Functions.calc_probe(rsample[i],20,R,Rbar,kg,Bi,Henry,DCO2,Cinit,PCO2)
# C_all_t30[i]=DropletDiffusion_Functions.calc_probe(rsample[i],30,R,Rbar,kg,Bi,Henry,DCO2,Cinit,PCO2)
# plt.figure()
# plt.plot(tsample,C_all_Rover4,label='R/4')
# plt.plot(tsample,C_all_Rover2,label='R/2')
# plt.plot(tsample,C_all_Rover1,label='R')
# plt.xlabel('time (s)')
# plt.ylabel('Concentration (mol/m3)')
# plt.legend(loc='best')
# plt.title('Time history of concentration at various locations of the droplet')
# # plt.savefig('%s/Concentration_timehistory.pdf'%WritetoFolder, bbox_inches='tight')
# plt.figure()
# plt.plot(rsample,C_all_t10,label='t=10s')
# plt.plot(rsample,C_all_t20,label='t=20s')
# plt.plot(rsample,C_all_t30,label='t=30s')
# plt.xlabel('r (m)')
# plt.ylabel('Concentration (mol/m3)')
# plt.legend(loc='best')
# plt.title('Radial distribution of concentration in the droplet at various times')
# # plt.savefig('%s/Concentration_timehistory.pdf'%WritetoFolder, bbox_inches='tight')
```
```python
# -*- coding: utf-8 -*-
"""
Created on Fri Feb 4 17:50:33 2022
@author: ab65462
"""
#### Includes design parameters, droplet size, add more..
import os
import math
import numpy as np
import pandas as pd
import scipy.interpolate
### Universal Constants-
g = 9.8 # Gravitational acceleration [m/s2]
Rbar = 8.3145 # Universal gas constant [J/mol/K]
### Pressure and Temperature at Inlet and ambient/operating conditions
P_tot = 101.325 #Ambient Pressure [kPa]
T = 25 + 273.15 #Inlet Temperature [K]
P_CO2_in = 15.000 #Inlet CO2 partial pressure [kPa]
X1= P_CO2_in/P_tot
Y1 = X1/(1-X1)
C_init = 5 #Initial Concentration in [M] or [kmol/m3]
###Reactor/Channel Design Specs-
D = 0.1 #Inner Diameter of Channel [m]
H = 0.55 #Height of channel [m]
A_D = np.pi*D*D/4 #Cross sectional area of channel
#nozzle_type = 'P-40'
#nozzle_type = 'P-28'
#nozzle_type = 'MPL30N'
folder = r"C:\Users\ab65462\Desktop\CO2 Musk Foundation\Spray Model\Kuntz_Aroonwilas_Data"
#fig_validate = '2b-1.91.csv'
#fig_validate = '2b-4.59.csv'
#fig_validate = '2fc-P28.csv'
fig_validate = '2f-P40.csv'
filename = os.path.join(folder,fig_validate)
df = pd.read_csv(filename)
kG_v_Q = df.to_numpy()
x_exp = np.array(kG_v_Q[:,0])
y_exp = np.array(kG_v_Q[:,1])
### Need to create filename for df so that these if conditions become cleaner
if fig_validate == '2b-1.91.csv':
nozzle_type = 'P-28'
Q_l_flux = 1.91 # Flux of liquid flow rate in m3/m2-h
Var = 'Gas Flow Rate'
if fig_validate == '2b-4.59.csv':
nozzle_type = 'P-28'
Q_l_flux = 4.59 # Flux of liquid flow rate in m3/m2-h
Var = 'Gas Flow Rate'
if fig_validate == '2fc-P28.csv':
nozzle_type = 'P-28'
Q_g_flux = 382 # Flux of gas flow rate in m3/m2-h
Var = 'Liquid Flow Rate'
if fig_validate == '2f-P40.csv':
nozzle_type = 'P-40'
Q_g_flux = 382 # Flux of gas flow rate in m3/m2-h
Var = 'Liquid Flow Rate'
#nozzle_type = 'MPL30N'
####Nozzle Specs-
if nozzle_type == 'P-40':
df_Q_vs_delP = pd.read_csv(r'C:\Users\ab65462\Desktop\CO2 Musk Foundation\Spray Model\Manufacturer Data\P-40-gpm_vs_psi.csv')
df_d_vs_delP = pd.read_csv(r'C:\Users\ab65462\Desktop\CO2 Musk Foundation\Spray Model\Manufacturer Data\P-40-micron_vs_psi.csv')
QvP = df_Q_vs_delP.to_numpy()
dvP = df_d_vs_delP.to_numpy()
d_orifice = 0.042*0.254 #[m]
K_factor = .0443 #Manufacturer's design parameter for nozzle that relates deltaP and GPM
if nozzle_type == 'P-28':
df_Q_vs_delP = pd.read_csv(r'C:\Users\ab65462\Desktop\CO2 Musk Foundation\Spray Model\Manufacturer Data\P-28-gpm_vs_psi.csv')
df_d_vs_delP = pd.read_csv(r'C:\Users\ab65462\Desktop\CO2 Musk Foundation\Spray Model\Manufacturer Data\P-28-micron_vs_psi.csv')
QvP = df_Q_vs_delP.to_numpy()
dvP = df_d_vs_delP.to_numpy()
d_orifice = 0.028*0.254 #[m]
K_factor = .0206 #Manufacturer's design parameter for nozzle that relates deltaP and GPM
if nozzle_type == 'MPL30N':
df_Q_vs_delP = pd.read_csv(r'C:\Users\ab65462\Desktop\CO2 Musk Foundation\Spray Model\Manufacturer Data\MPL30N-gpm_vs_psi.csv')
df_d_vs_delP = pd.read_csv(r'C:\Users\ab65462\Desktop\CO2 Musk Foundation\Spray Model\Manufacturer Data\MPL30N-micron_vs_psi.csv')
QvP = df_Q_vs_delP.to_numpy()
dvP = df_d_vs_delP.to_numpy()
d_orifice = 1.1*1e-3 #[m]
#K_factor = ? #Manufacturer's design parameter for nozzle that relates deltaP and GPM
def d_from_Q(Q_GPM): #Returns droplet diameter in um from flowrate in gpm
f_PvQ = scipy.interpolate.interp1d(QvP[:,1], QvP[:,0])
f_dvP = scipy.interpolate.interp1d(dvP[:,0], dvP[:,1])
delP = f_PvQ(Q_GPM)
#print(delP, 'psi')
droplet_dia = f_dvP(delP)
return droplet_dia
#print(d_from_Q(.3))
```
```python
# -*- coding: utf-8 -*-
"""
Created on Wed Apr 6 21:05:44 2022
@author: ab65462
"""
### Includes simple side functions like viscosity, henry constant etc calculators
import numpy as np
import math
import matplotlib.pyplot as plt
import scipy.interpolate
from Input_Params import *
from scipy.integrate import *
def D_and_H_CO2_l_calc(T): # Temp input is in K, returns D and H for CO2 in 5M MEA at 30C
D_N2O_w = 5.07*1e-6*math.exp(-2371/T) # m2/s ADDD REFERENCE FOR THESE if it is not Tsai etal 2000!!!!
D_CO2_w = 2.35*1e-6*math.exp(-2119/T) # m2/s
D_N2O_MEA = 1.41*1e-9 # 30C and 5M MEA, data from Jium-Jie Ko et al 2001, [m2/s]
D_CO2_MEA = D_N2O_MEA*D_CO2_w/D_N2O_w
H_N2O_w = 8.547*1e6*math.exp(-2284/T) #[kPa-m3/kmol]
H_CO2_w = 2.8249*1e6*math.exp(-2044/T) #[kPa-m3/kmol]
H_N2O_MEA = 4924 #30C and 5M, data from Tsai et al 2000, [kPa-m3/kmol]
H_CO2_MEA = H_N2O_MEA*H_CO2_w/H_N2O_w
##### need to add a function of temp to H_N2O_MEA and D_N2O_MEA, then later add effect of concentration on all this too
return D_CO2_MEA, H_CO2_MEA
def k2_calc(T): #Temp input is in K
k2 = 1.72*1e11*math.exp(-4915.5/T) #Ramezaini etal 2021, 2nd order reaction rate used here, [m3/mols]
return k2
def calc_viscosity(T): #T in Kelvin, and dynamic viscosity in [kg/m-s]
#The following relation is based on Sutherland's formula for 3 coefficients (values are taken for air) ...ref- https://www.afs.enea.it/project/neptunius/docs/fluent/html/ug/node294.htm
T0 = 273.11 #K
mu_0 = 1.716*1e-5 #kg/m-s
Suth = 110.56 #Sutherland Constant (K)
mu = mu_0*((T/T0)**1.5)*((T0+Suth)/(T+Suth))
return mu
#plt.figure()
def U_calc(d, Ujet, nu_gas, U_G, H, rho_gas, rho_drop):
tau = (d**2)*(rho_drop/rho_gas)/(18*nu_gas)
#print(tau)
#print((U_G+Ujet)*d/nu_gas, 'Re')
U0 = Ujet + U_G # Adjusts the jet velocity to relative velocity wrt air (counterflow)
n_points = 1001
tf = 5
t = np.linspace(0,tf,n_points)
U = np.zeros(n_points) # Absolute Velocity of fluid in m/s
s = np.zeros(n_points)
reached_column_end = 0
for i in range(n_points):
U[i] = U0*math.exp(-t[i]/tau) + g*tau*(1-math.exp(-t[i]/tau)) - U_G
s[i] = U0*tau*(1-math.exp(-t[i]/tau)) + g*tau*t[i] - g*(tau**2)*(1-math.exp(-t[i]/tau)) - U_G*t[i]
#print(t[i], s[i], U[i])
if s[i] > H and reached_column_end == 0:
print('Here!!')
reached_column_end = 1
n_points_needed = i+1
residence_time = (t[i] +t[i-1])/2
plt.plot(t[:], s[:], 'b-', label = 'U vs s')
#print(s[i])
#residence_time = (t[i-1])/2
print(residence_time, 't_res')
break
n_points_needed = i+1
print(n_points_needed)
plt.plot(t[:], s[:], 'g-', label = 'U vs s')
U_fn_t = np.zeros(n_points_needed)
s_fn_t = np.zeros(n_points_needed)
t_fn_t = np.zeros(n_points_needed)
for j in range(n_points_needed):
U_fn_t[j] = U[j]
s_fn_t[j] = s[j]
t_fn_t[j] = t[j]
z = np.linspace(0,H,n_points_needed)
U_z = scipy.interpolate.interp1d(s_fn_t, U_fn_t)
U_fn_z = U_z(z)
t_z = scipy.interpolate.interp1d(s_fn_t, t_fn_t)
t_fn_z = t_z(z)
#s_z = scipy.interpolate.interp1d(s_fn_t, U_fn_t)
#plt.plot(z[:], U_fn_z[:], 'go', label = 'U vs s')
return U_fn_z, s_fn_t, t_fn_z, residence_time, z
def calc_gas_prop(P_CO2_in,T):
X1 = P_CO2_in/P_tot
mu_air=calc_viscosity(T) # named as mu_air because it is currently evaluating for air only and not for flue gas
mu_gas = mu_air # Assumed viscosity of flue gas is same as that of air
M_gas = 0.028*(1-X1) + 0.044*(X1)
rho_gas = P_tot*1000/(Rbar/M_gas)/(T)
nu_gas = mu_gas/rho_gas
return nu_gas, rho_gas
def U_calc2(Q_l_flux, Q_g_flux, T):
Q_L = Q_l_flux*A_D/3600 #Liquid flow rate [m3/s] when Q_l_flux is in m3/m2-h
Q_L_GPM = Q_L*15850.3 #Liquid Flow Rate in Gallons/minute
d = d_from_Q(Q_L_GPM)*1e-6 #droplet SMD in [m]
Q_G =Q_g_flux*A_D/3600 # Gas Flow Rate [m3/s] (For aroonwilas paper, the flow rate is assumed total gas flow rate for now as info is not clerarly mentioned in paper)
U_G = Q_G/A_D # Gas Velocity [m/s]
U_jet = Q_L*4/(3.14*d_orifice*d_orifice) # Velocity of the droplet.. jet velocity is used here [m/s]
nu_gas, rho_gas = calc_gas_prop(P_CO2_in,T)
rho_drop = 1000 # Density of the droplet [kg/m3]. Both MEA and water have a density around 1000
U_fn_z, s_fn_z, t_fn_z, t_res, z = U_calc(d, U_jet, nu_gas, U_G, H, rho_gas, rho_drop)
U = H/t_res #Average velocity
return U
def hm_calc(Q_l_flux, Q_g_flux, P_CO2_in, T, CO2_loading, U):
Q_L = Q_l_flux*A_D/3600 #Liquid flow rate [m3/s] when Q_l_flux is in m3/m2-h
Q_L_GPM = Q_L*15850.3 #Liquid Flow Rate in Gallons/minute
d = d_from_Q(Q_L_GPM)*1e-6 #droplet SMD in [m]
R = d/2
Vol_1d = (4/3)*3.14*(d**3)/8
X1 = P_CO2_in/P_tot # Mole fraction of CO2 at inlet
####some gas quantities are for air, rest for flue gas, I need to incorporate the air part to flue gas as well later...
mu_air=calc_viscosity(T) # named as mu_air because it is currently evaluating for air only and not for flue gas
mu_gas = mu_air # Assumed viscosity of flue gas is same as that of air
M_gas = 0.028*(1-X1) + 0.044*(X1)
rho_gas = P_tot*1000/(Rbar/M_gas)/(T)
nu_gas = mu_gas/rho_gas
nu_gas, rho_gas = calc_gas_prop(P_CO2_in,T)
######remove####DCO2_w=1.92e-9 # Diffusivity of CO2 in the droplet [m2/s]. This value is for MEA wt30%
D_CO2_air = 1.6e-5 # Diffusivity of CO2 in the air [m2/s]
######remove#####Henry=4219 # Henry's coefficient in [m3Pa/mol]. This value is for water (AbuArabi 2001)
rho_drop = 1000 # Density of the droplet [kg/m3]. Both MEA and water have a density around 1000
mu_drop = 2.48*1e-3 # Dynamic viscosity of the droplet at about 30C and 30wt% or 5M MEA [Pas] Ref: Arachchigeetal 2013
nu_drop = mu_drop/rho_drop
D_CO2_MEA, H_CO2_MEA = D_and_H_CO2_l_calc(T)
#print(D_CO2_MEA, H_CO2_MEA)
#Cinit=0 # Initial CO2 concentration in the droplet [mol/m3] .. To be used later when C(z) will be added
k2 = k2_calc(T) #m3/mol-s
C_drop = C_init #kmol/m3 or M
#print(k2)
print(CO2_loading)
1/0
# if CO2_loading > 0.5:
# if CO2_loading < 0.5:
k_liq_w_kin = ((k2*D_CO2_MEA*C_init*(1-2*CO2_loading)*1000)**0.5/H_CO2_MEA) #[mol/(Pa-m2-s)] or [mol-s/(kg-m)]
#print(k_liq_w_kin)
###Flow Parameters-
# Q_G =Q_g_flux*A_D/3600 #Gas Flow Rate [m3/s] (For aroonwilas paper, the flow rate is assumed total gas flow rate for now as info is not clerarly mentioned in paper)
# U_G = Q_G/A_D #Gas Velocity [m/s]
U_jet = Q_L*4/(3.14*d_orifice*d_orifice) # Velocity of the droplet.. jet velocity is used here [m/s]
# # U_fn_z, s_fn_z, t_fn_z, t_res, z = U_calc(d, U_jet, nu_gas, U_G, H, rho_gas, rho_drop) ### Gives velocity of droplet at each location of height, and droplet residence time
# U = H/t_res #Average velocity
# #U = U_jet
#print('Gas velocity is', U_G, 'and droplet velocity is', U)
# tau = d*d/(18*nu_drop) #Time constant for velocity variation
# U_rel_terminal = g*tau
# U_terminal = U_rel_terminal-U_G
# #print( 'Gas velocity is', U_G, 'and droplet velocity is', U, 'and terminal velocity is', U_terminal)
# U_array[j] = U_terminal
# print(U_terminal, U, U_G)
Re_drop = U*d/nu_gas ###Should it be relative velocity or absolute fluid velocity
Sc_drop = nu_gas/D_CO2_air # Check these!! Also check if the diffusion represented here is co2 or amine
Sh_drop = 2 + 0.552*(Re_drop**(1/2))*(Sc_drop**(1/3)) #Froessling Equation
#print(Sh_drop, 'is sherwood number')
#print('nu_gas=',nu_gas)
Re_jet = U_jet*d/nu_gas
Stop_dist = (rho_drop*d/rho_gas)*(Re_jet**(1/3)-math.sqrt(6)*math.atan(Re_jet**(1/3)/math.sqrt(6)))
#print('S=', Stop_dist)
kg=Sh_drop*D_CO2_air/(Rbar*T*d) # [mol/(Pa-m2-s)] or [mol-s/(kg-m)]
#print(kg)
#kG = kg
kG = 1/(1/kg + 1/k_liq_w_kin)
hm = kG*Rbar*T #[m/s]
# print('kL = ', k_liq_w_kin, 'kg = ', kg, 'and kG = ', kG)
#print('Kgav in mol/Pa-m3-hr assuming area for all droplets = ', kG*3600*H*Q_L/(U*Vol_1d*A_D*H) )
#print(hm, 'is hm')
Bi=kg*(d)*Rbar*T/D_CO2_air
#print("Biot is %2.2f" %(Bi))
return hm, d
```
```
```
Sign in with Wallet
Connect another wallet