James - Script

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])
################################################################

#BED
NUMBER_BEDS = 10 #Physical packed beds
DIAMETER_BED = 0.0779 #metres
HEIGHT_BED = 0.27 #metres
AREA_BED = ((math.pi)/4)*DIAMETER_BED**2 #m^2
VOLUME_BED = AREA_BED*HEIGHT_BED #m^3
V_TOTAL = VOLUME_BED * NUMBER_BEDS #total beds volume (m^3)

#PARTICLE
DIAMETER_PARTICLE = 0.005 #metres
VOLUME_PARTICLE = (1/6)*math.pi*(DIAMETER_PARTICLE**3) #particle volume m^3
PARTICLE_NUMBER = math.floor(V_TOTAL/VOLUME_PARTICLE) #total number of particles
SA_PARTICLE = 4*math.pi*((DIAMETER_PARTICLE/2)**2) #m^2
SA_TOTAL = SA_PARTICLE*PARTICLE_NUMBER #total surface area (m^2)
RHO_ALUMINA = 3950 #kg/m3
M_ALUMINA_KG = VOLUME_PARTICLE*PARTICLE_NUMBER*RHO_ALUMINA # total mass of alumina (kg)
MW_ALUMINA = 101.96 #g/mol
M_ALUMINA_MOL = M_ALUMINA_KG*1000/MW_ALUMINA #total mass of alumina (mol)
CP_ALUMINA = 89.72 #J/mol.K
K_ALUMINA = 5 #W/m.K

#FLUID
MW_AIR = 29 #g/mol

#INITIAL CONDITIONS
T0 = 80  # initial fluid temperature (K)
Ts0 = 80  # initial solid temperature (K)
P = 3000000  # Pa
H0 = flasher.flash(T=T0, P=P, zs=zs, retry = True).H() # bed fluid starting enthalpy (J/mol)


#BOUNDARY CONDITIONS
Tin = 300  # Inlet fluid temperature (K)
minkg = 0.022 #kg/s
min = minkg*1000/MW_AIR #mol/s
hIn = flasher.flash(T=Tin, P=P, zs=zs, retry = True).H()  # Inlet fluid enthalpy (J/mol)

#SPATIAL DEPENDENCE OF QUANTITIES
n = 3 # number of elements in spatial array
V = V_TOTAL/n
sa = SA_TOTAL/n #surface area per unit volume
mAlumina = M_ALUMINA_MOL/n #alumina mass per unit volume

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

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

def dpdh(H, P):
    #returns the 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_dot returns U W/m2.K
    air = Mixture(['nitrogen', 'oxygen'], Vfgs=[0.79, 0.21], H=H, P=P)
    mu = air.mu*1000/MW_AIR #mol/m.s
    cp = air.Cpm #J/mol.K
    kg = air.k #W/m.K
    g0 = m/AREA_BED #mol/m2.s
    a = sa*n/V_TOTAL #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)/DIAMETER_PARTICLE

    h = 1/((1/hfs) + ((DIAMETER_PARTICLE/beta)/K_ALUMINA))

    return h

# def vapourfraction(H,P):
#     x = flasher.flash(H=H, P=P, zs=zs, retry=True).VF # mol basis
#     return x

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


def my_system(t, y, n, hIn, min, mAlumina, CP_ALUMINA, sa, V):
    dydt = np.zeros(2 * n) #setting up zeros array

    i=0

    ## Using thermo
    T = temperature(y[i],P) #initial fluid temperature
    m = mass(y[i],P) #initial mass holdup
    print("T = ", T)
    print("M = ", m)

    #initial boundary values
    dydt[i] = (min * (hIn - y[i]) + (U(y[i],P,min) * sa * (y[i + 1] - T))) / m # dH/dt (eq. 2)
    dydt[i + 1] = -(U(y[i],P,min) * sa * (y[i + 1] - T)) / (mAlumina * CP_ALUMINA) # dTs/dt from eq.3
    dmdt = dydt[i] * dpdh(y[i], P) * V # dm/dt eq. 4b
    mj = min - dmdt # mass flow out (eq.4a)

    #Solution array y will look like [h0,Ts0,mdot0, h1, Ts1, mdot1,...hn,Tsn,mdotn] at all values of t
    for i in range(2, 2 * n, 2): #starting at index 2, and incrementing by 2 because we are solving for 'triplets' [h,Ts] in each loop

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

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

    return dydt


y0 = np.copy(np.asarray([H0, Ts0] * n)) #setting up initial condition array (H0,Ts0 and m at all points
y0[0] = H0 # initial fluid bed enthalpy
y0[1] = Ts0 # initial solid bed temperature


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

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

h_arr = solution.y[0::2]
Ts_arr = solution.y[1::2]

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

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

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

dmdt_arr = np.zeros_like(dhdt_arr) #change in mass holdup
T_arr = np.zeros_like(dhdt_arr) # fluid temperature
m_hold_arr = np.zeros_like(dhdt_arr) # mass hold up
U_arr = np.zeros_like(dhdt_arr) # heat transfer coefficient array
mj_arr = np.zeros_like(dhdt_arr) # heat transfer coefficient array
x_arr = np.zeros_like(dhdt_arr)

for i in range(time.shape[0]): #iteratung through 2D array time x #positions
    for j in range(n): #yep

        T = temperature(h_arr[j,i],P)

        dmdt_arr[j, i] = dhdt_arr[j, i] * dpdh(h_arr[j, i], P) * V
        T_arr[j,i] = T
        m_hold_arr[j, i] = mass(h_arr[j, i], P)
        U_arr[j,i] = mAlumina*CP_ALUMINA*dTsdt_arr[j,i]/(sa*(T_arr[j,i] - Ts_arr[j,i]))
        #x_arr[j,i] = vapourfraction(h_arr[j,i], P)


# creating Tf_arr from Ts_arr
Tf_arr = np.zeros_like(dhdt_arr)
Tf_arr = Ts_arr + dTsdt_arr * (mAlumina * CP_ALUMINA) / (U_arr * sa)
# Tf_arr[:, 0] = Ts0  # I.C cannot be calculated as there is no derivative at the boundary. Manually entering here

print("simulation complete")

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

plt.ylabel(r'Fluid temperature')
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 (mol)')
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,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)')
plt.xlabel(r'Time')
# plt.ylim((80, 91))
plt.grid('both')
plt.title("h for all tanks")
plt.legend(loc='lower right')

plt.show()