fluid_packed_bed

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from chemicals import normalize
from thermo import ChemicalConstantsPackage, CEOSGas, CEOSLiquid, PRMIXTranslatedConsistent, FlashVL
from thermo.interaction_parameters import IPDB
import math
from thermo.chemical import Mixture

##Setup of the fluid mixture properties (air)
constants, properties = ChemicalConstantsPackage.from_IDs(['nitrogen', 'oxygen', 'argon'])
kijs = IPDB.get_ip_asymmetric_matrix('ChemSep PR', constants.CASs, 'kij')
eos_kwargs = {'Pcs': constants.Pcs, 'Tcs': constants.Tcs, 'omegas': constants.omegas, 'kijs': kijs}
gas = CEOSGas(PRMIXTranslatedConsistent, eos_kwargs=eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases)
liquid = CEOSLiquid(PRMIXTranslatedConsistent, eos_kwargs=eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases)
flasher = FlashVL(constants, properties, liquid=liquid, gas=gas)
zs = normalize([78.08, 20.95, .93])
################################################################

numberBeds = 10
diameterBed = 0.0779 #metres
heightBed = 0.27
areaBed = ((math.pi)/4)*diameterBed**2
volumeBed = areaBed*heightBed
vTotal = volumeBed * numberBeds #total bed volume (m^3)

diameterParticle = 0.005
volumeParticle = (1/6)*math.pi*(diameterParticle**3)
particleNumber = math.floor(vTotal/volumeParticle)
saParticle = 4*math.pi*((diameterParticle/2)**2)
saTotal = saParticle*particleNumber #total surface area (m^2)

rhoAlumina = 3950
mAluminakg = volumeParticle*particleNumber*rhoAlumina # total mass of alumina (kg)
mWAlumina = 101.96
mAluminaTotal = mAluminakg*1000/mWAlumina #total mass of alumina (mol)

cpAlumina = 89.72
kAlumina = 5

mWAir = 29

minkg = 0.0022*10
min = minkg*1000/mWAir #mol/s

T0 = 80  # initial fluid temp (K)
Ts0 = 80  # initial solid temp (K)

Tin = 300  # Inlet fluid temperature (K)
P = 3000000  # Pa

n = 3 # number of elements in spatial array
V = vTotal/n
sa = saTotal/n #surface area per unit volume
mAlumina = mAluminaTotal/n #alumina mass per unit volume


## Using thermo
hIn = flasher.flash(T=Tin, P=P, zs=zs, retry = True).H()  # Inlet fluid enthalpy (J/mol)
H0 = flasher.flash(T=T0, P=P, zs=zs, retry = True).H() # bed fluid starting enthalpy (J/mol)

def temp(H,P):
    # returns temperature given enthalpy (after processing function)
    T = flasher.flash(H=H, P=P, zs=zs, retry=True).T
    return T

def mass(H, P):
    # returns mass holdup in mol
    m = flasher.flash(H=H, P=P, zs=zs, retry=True).rho()*V
    return m

def dpdh(H, P):
    #derivative of density w.r.t enthalpy at constant pressure
    res = flasher.flash(H=H, P=P, zs=zs, retry=True)
    if res.phase_count == 1:
        if res.phase == 'L':
            drho_dTf = res.liquid0.drho_dT()
        else:
            drho_dTf = res.gas.drho_dT()
    else:
        drho_dTf = res.bulk._equilibrium_derivative(of='rho', wrt='T', const='P')
    dpdh = drho_dTf/res.dH_dT_P()
    return dpdh

def U(H,P,m):
    # Given T, P, m, calculate heat transfer coefficient (function of T,P, mass flow)
    air = Mixture(['nitrogen', 'oxygen'], Vfgs=[0.79, 0.21], H=H, P=P)
    mu = air.mu*1000/mWAir #mol/m.s
    cp = air.Cpm #J/mol.K
    kg = air.k #W/m.K
    g0 = m/areaBed #mol/m2.s
    a = sa*n/vTotal #m^2/m^3 #QUESTIONABLE
    psi = 1
    beta = 10
    pr = (mu*cp)/kg
    re = (6*g0)/(a*mu*psi)
    hfs = ((2.19*(re**1/3)) + (0.78*(re**0.619)))*(pr**1/3)*(kg)/diameterParticle

    h = 1/((1/hfs) + ((diameterParticle/beta)/kAlumina))

    return h

# time steps
t_max = 100
t = (0, t_max)  # time range
# t_eval = np.linspace(0, t_max, 1000 * t_max)  # optional


def my_system(t, y, n, hIn, min, mAlumina, cpAlumina, sa, V):

    #setup of differential equations in algebraic form

    dydt = np.zeros(3 * n) #setting up zeros array for solution (solving for [H0,Ts0,m0,H1,Ts1,m1,H2,Ts2,m2,..Hn,Tsn,mn])

    # y = [h_0, Ts_0, m_0, ... h_n, Ts_n, m_n]
    # y[0] = hin
    # y[1] = Ts0
    # y[2] = minL

    i=0

    ## Using thermo
    T = temp(y[i],P) #initial T
    m = mass(y[i],P) #initial m

    #initial values
    dydt[i] = (min * (hIn - y[i]) + (U(hIn,P,min) * sa * (y[i + 1] - T))) / m # dH/dt (eq. 2)
    dydt[i + 1] = -(U(hIn,P,min) * sa * (y[i + 1] - T)) / (mAlumina * cpAlumina) # dTs/dt from eq.3
    dmdt = dydt[i] * dpdh(y[i], P) * V # dm/dt (holdup variation) eq. 4b
    dydt[i + 2] = min - dmdt # mass flow out (eq.4a)
    # print("H fluid 0: ", dydt[i])
    # print("Ts 0: ", dydt[i+1])
    # print("m dot: ", dydt[i+2])


    for i in range(3, 3 * n, 3): #starting at index 3, and incrementing by 3 because we are solving for 'triplets' [h,Ts,m] in each loop

        ## Using thermo
        T = temp(y[i],P)
        m = mass(y[i],P)
        # print("T check: ", T)

        # [h, TS, mdot]
        dydt[i] = (dydt[i-1] * (y[i - 3] - y[i]) + (U(y[i-3], P, dydt[i-1]) * sa * (y[i + 1] - T))) /  m # dH/dt (eq.2), dydt[i-1] is the mass of the previous tank
        dydt[i + 1] = -(U(y[i-3], P, dydt[i-1]) * sa * (y[i + 1] - T)) / (mAlumina * cpAlumina) # dTs/dt eq. (3)
        dmdt = dydt[i] * dpdh(y[i], P) * V # Equation 4b
        dydt[i + 2] = dydt[i-1] - dmdt # Equation 4a

        # print("H fluid: ", dydt[i])
        # print("Ts: ", dydt[i + 1])
        # print("m dot: ", dydt[i + 2])
    return dydt

y0 = np.copy(np.asarray([H0, Ts0, min] * n)) #setting up initial condition array (H0,Ts0 and m_dot at all points
print("yo")
print(y0)
y0[0] = H0 # initial fluid bed enthalpy
y0[1] = Ts0 # initial solid bed temperature
y0[2] = min #initial mass holdup

parameters = (n, hIn, min, mAlumina, cpAlumina, sa, V)

solution = solve_ivp(my_system, t, y0, args=parameters, method='Radau')

##FROM HERE TO BELOW IS AFTER PROCESSING (EXTRACTION OF SOLUTION ARRAYS, PLOTTING, ETC)
h_arr = solution.y[0::3]
Ts_arr = solution.y[1::3]
m_arr = solution.y[2::3]

time = solution.t
# dt = 2*(time[2:]-time[:-2])

dydt = np.zeros((3 * n, time.size))
for i in range(time.size):
    dydt[:, i] = my_system(t, solution.y[:, i], n, hIn, min, mAlumina, cpAlumina, sa, V)

dhdt_arr = dydt[0::3]
dTsdt_arr = dydt[1::3]
mdot_arr = dydt[2::3]

dmdt_arr = np.zeros_like(mdot_arr) #
T_arr = np.zeros_like(mdot_arr) # fluid temperature
m_hold_arr = np.zeros_like(mdot_arr) # mass hold up
U_arr = np.zeros_like(mdot_arr) # mass hold up

for i in range(time.shape[0]):
    for j in range(n):
        # T = temp(h_arr[j,i])

        ## Using thermo
        T = temp(h_arr[j,i],P)

        T_arr[j,i] = T

        ## Using thermo
        m_hold_arr[j,i] = mass(h_arr[j,i],P)

        dmdt_arr[j,i] = dhdt_arr[j,i] * dpdh(h_arr[j,i], P) * V

        U_arr[j,i] = mAlumina*cpAlumina*dTsdt_arr[j,i]/(sa*(T_arr[j,i] - Ts_arr[j,i]))

# creating the fluid array from the solid array
Tf_arr = np.zeros_like(Ts_arr)
Tf_arr = Ts_arr + dTsdt_arr * (mAlumina * cpAlumina) / (U_arr * sa)
# Tf_arr[:, 0] = Ts0  # I.C cannot be calculated as there is no derivative at the boundary. Manually entering here
#data1 = np.concatenate([np.array([time]), h_arr]).transpose()
# data1 = np.stack((data1,data1),axis=2)

#np.savetxt("data1.txt",data1)
print("simulation omplete")


#print(x_arr)


'''
snapshot = 1000  # snapshot time spacing: 1000=1s
style = ['-', '--', '-.', ':', '-','-', '--', '-.', ':', '-','-', '--']
color = ['b','r','g','c','b','r','g','c','b','r','g','c']
for i in range(len(style)):
    plt.plot(mdot_arr[:, i*snapshot], style[i],
             color=color[i],
             label="t="+str(i*snapshot/1000)+"s")

plt.ylabel(r'm (kg/s)')
plt.xlabel(r'position(j)')
plt.legend(loc='lower left')
plt.show()
'''
# plt.plot(time,h_arr[0, :], '-.', color='b', label="position 0 h")
# plt.plot(time,h_arr[1, :], '--', color='b', label="position 1")
# plt.plot(time,h_arr[2, :], ':', color='b', label="position 2")
# plt.plot(time,Tf_arr[0, :], '-.', color='g', label="position 0 T")
#plt.plot(time,x_arr[0, :], '-.', color='r', label="tank 1 quality")
#plt.plot(time,x_arr[1, :], '-.', color='g', label="tank 2 quality")
#plt.plot(time,x_arr[2, :], '-.', color='b', label="tank 3 quality")
'''
plt.figure()
plt.plot(time,dmdt_arr[0, :], '-', color='r', label="tank 1 dm/dt")
plt.plot(time,dmdt_arr[1, :], '--', color='g', label="tank 2 dm/dt")
plt.plot(time,dmdt_arr[2, :], '-.', color='b', label="tank 3 dm/dt")
plt.plot(time,dmdt_arr[3, :], '--', color='m', label="tank 4 dm/dt")
plt.plot(time,dmdt_arr[4, :], '-.', color='c', label="tank 5 dm/dt")

# plt.plot(time,Tf_arr[1, :], '--', color='g', label="position 1")
# plt.plot(time,Tf_arr[2, :], ':', color='g', label="position 2")
# from plot_together import plot_together
# plot_together(data1[:,[0,1]:], '--', color='b', label="position 0")
# plot_together(data1[:,[0,2],:], '--', color='b', label="position 1")
# plot_together(data1[:,[0,3],:], ':', color='b', label="position 2")
plt.ylabel(r'dm/dt')
plt.xlabel(r'Time')
# plt.ylim((80, 91))
plt.grid('both')
plt.title("dm/dt for all tanks")
plt.legend(loc='lower right')
'''

plt.figure()
plt.plot(time,mdot_arr[0, :], '-', color='b', label="tank 1 mdot")
plt.plot(time,mdot_arr[1, :], '--', color='b', label="tank 2 mdot")
plt.plot(time,mdot_arr[2, :], '-.', color='b', label="tank 3 mdot")
#plt.plot(time,mdot_arr[3, :], '--', color='b', label="tank 4 mdot")
#plt.plot(time,mdot_arr[4, :], '-.', color='b', label="tank 5 mdot")

plt.ylabel(r'Mass Flow (kg/s)')
plt.xlabel(r'Time')
# plt.ylim((80, 91))
plt.grid('both')
plt.title("m dot for all tanks")
plt.legend(loc='lower right')

#########################################################
plt.figure()
plt.plot(time,m_arr[0, :], '-', color='b', label="tank 1  m_flow")
plt.plot(time,m_arr[1, :], '--', color='b', label="tank 2 m_flow")
plt.plot(time,m_arr[2, :], '-.', color='b', label="tank 3 m_flow")
#plt.plot(time,m_arr[3, :], '--', color='b', label="tank 4 m_flow")
#plt.plot(time,m_arr[4, :], '-.', color='b', label="tank 5 m_flow")

plt.ylabel(r'Mass flow')
plt.xlabel(r'Time')
# plt.ylim((80, 91))
plt.grid('both')
plt.title("Mass flow for all tanks")
plt.legend(loc='lower right')

###########################################################
plt.figure()
plt.plot(time,T_arr[0, :], '-', color='b' , label="tank 1 Fluid temp", marker = "x")
plt.plot(time,T_arr[1, :], '--', color='b', label="tank 2 Fluid temp")
plt.plot(time,T_arr[2, :], '-.', color='b', label="tank 3 Fluid temp")
#plt.plot(time,T_arr[3, :], '--', color='b', label="tank 4 Fluid temp")
#plt.plot(time,T_arr[4, :], '-.', color='b', label="tank 5 Fluid temp")

plt.ylabel(r'Fluid temp')
plt.xlabel(r'Time')
# plt.ylim((80, 91))
plt.grid('both')
plt.title("Fluid temperature for all tanks")
plt.legend(loc='lower right')

###########################################################
plt.figure()
plt.plot(time,m_hold_arr[0, :], '-', color='b' , label="tank 1 Mass hold up")
plt.plot(time,m_hold_arr[1, :], '--', color='b', label="tank 2 Mass hold up")
plt.plot(time,m_hold_arr[2, :], '-.', color='b', label="tank 3 Mass hold up")
#plt.plot(time,m_hold_arr[3, :], '--', color='b', label="tank 4 Mass hold up")
#plt.plot(time,m_hold_arr[4, :], '-.', color='b', label="tank 5 Mass hold up")

plt.ylabel(r'Mass hold up')
plt.xlabel(r'Time')
# plt.ylim((80, 91))
plt.grid('both')
plt.title("Mass hold up for all tanks")
plt.legend(loc='lower right')

###########################################################
plt.figure()
plt.plot(time,mdot_arr[0, :], '-', color='b' , label="tank 1 Mass flow rate")
plt.plot(time,mdot_arr[1, :], '--', color='b', label="tank 2 Mass flow rate")
plt.plot(time,mdot_arr[2, :], '-.', color='b', label="tank 3 Mass flow rate")
#plt.plot(time,mdot_arr[3, :], '--', color='b', label="tank 4 Mass flow rate")
#plt.plot(time,mdot_arr[4, :], '-.', color='b', label="tank 5 Mass flow rate")

plt.ylabel(r'Mass flow rate')
plt.xlabel(r'Time')
# plt.ylim((80, 91))
plt.grid('both')
plt.title("Mass flow rate for all tanks")
plt.legend(loc='lower right')

###########################################################
plt.figure()
plt.plot(time,Ts_arr[0, :], '-', color='b' , label="tank 1 Solid Temperature")
plt.plot(time,Ts_arr[1, :], '--', color='b', label="tank 2 Solid Temperature")
plt.plot(time,Ts_arr[2, :], '-.', color='b', label="tank 3 Solid Temperature")
#plt.plot(time,Ts_arr[3, :], '--', color='b', label="tank 4 Solid Temperature")
#plt.plot(time,Ts_arr[4, :], '-.', color='b', label="tank 5 Solid Temperature")

plt.ylabel(r'Solid Temperature')
plt.xlabel(r'Time')
# plt.ylim((80, 91))
plt.grid('both')
plt.title("Solid Temperature for all tanks")
plt.legend(loc='lower right')

###########################################################
plt.figure()
plt.plot(time,h_arr[0, :], '-', color='b', label="tank 1 h")
plt.plot(time,h_arr[1, :], '--', color='b', label="tank 2 h")
plt.plot(time,h_arr[2, :], '-.', color='b', label="tank 3 h")
#plt.plot(time,h_arr[3, :], '--', color='b', label="tank 4 mdot")
#plt.plot(time,h_arr[4, :], '-.', color='b', label="tank 5 mdot")

plt.ylabel(r'h (J/mol.K)')
plt.xlabel(r'Time')
# plt.ylim((80, 91))
plt.grid('both')
plt.title("h for all tanks")
plt.legend(loc='lower right')

plt.show()