fluid model

import numpy as np
import scipy as sp
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import math
# from stream import Stream
from thermo.mixture import Mixture

################################################################
#IGNORE THIS BLOCK OF CODE IT IS JUST THERE TO SET UP A MIXTURE OBJECT (that allows thermodynamic property methods to be called)
################################################################
from chemicals import normalize
from thermo import ChemicalConstantsPackage, CEOSGas, CEOSLiquid, PRMIXTranslatedConsistent, FlashVL
from thermo.interaction_parameters import IPDB, SPDB
constants, properties = ChemicalConstantsPackage.from_IDs(['nitrogen', 'oxygen', 'argon'])
kijs = IPDB.get_ip_asymmetric_matrix('ChemSep PR', constants.CASs, 'kij')
Ls = SPDB.get_parameter_vector(name='PRTwu_PinaMartinez', CASs=constants.CASs, parameter='TwuPRL')
Ms = SPDB.get_parameter_vector(name='PRTwu_PinaMartinez', CASs=constants.CASs, parameter='TwuPRM')
Ns = SPDB.get_parameter_vector(name='PRTwu_PinaMartinez', CASs=constants.CASs, parameter='TwuPRN')
cs = SPDB.get_parameter_vector(name='PRTwu_PinaMartinez', CASs=constants.CASs, parameter='TwuPRc')
alpha_coeffs = [(Ls[i], Ms[i], Ns[i]) for i in range(constants.N)]
eos_kwargs = {'Pcs': constants.Pcs, 'Tcs': constants.Tcs, 'omegas': constants.omegas, 'kijs': kijs,
'cs': cs, 'alpha_coeffs':alpha_coeffs}
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])
################################################################

n = 3  # number of elements in spatial array

volumeV = 0.01 / n  # m3
area = 0.1 / n  # m2
mS = 1 / n  # kg

P = 200000  # Paß


mS = 57.75 / n  # kg
cpS = 1.1  # kJ/kgK
u = 1  # heat transfer coefficient between fluid and solid (W/m2.K)
T0 = 80  # initial fluid temp (K)

Tin = 300  # Inlet fluid temperature (K)
minL = 0.05  # inlet mass flow (kg/s)

## Using thermo
hIn = flasher.flash(T=Tin, P=P, zs=zs).H_mass()/1000  # Inlet enthalpy (J/kg)
H0 =  flasher.flash(T=T0, P=P, zs=zs).H_mass() /1000 # fluid starting enthalpy (J/kg)

print("thermo hIn = ", hIn)
print("thermo H0 = ", H0)

global equil_state

## First equilibrium state to accelerate convergence
equil_state = flasher.flash(H=H0, P=P, zs=zs, retry=True)

def temp_thermo(H,P,last_equil_state=None):
    # returns temperature given enthalpy (after processing function)
    equil_state = flasher.flash(H=H, P=P, zs=zs, retry=True,hot_start=last_equil_state)
    T_n = equil_state.T
    print("H = ",H)
    print("T = ", T_n)
    return T_n, equil_state


def mass_thermo(H, P, last_equil_state=None):
    # returns mass holdup given enthalpy
    equil_state = flasher.flash(H=H, P=P, zs=zs, retry=True,hot_start=last_equil_state)
    m = equil_state.rho()*volumeV
    return m, equil_state


def dpdh(H, P):

    ## Using thermo
    res = flasher.flash(H=H, P=P, zs=zs, retry=True)
    if res.phase_count == 1:
        if res.phase == 'L':
            print("liquid phase")
            drho_dT = 0
        else:
            drho_dT = res.gas.drho_dT_H()
            print("gas phase")
    else:
        drho_dT = res.bulk._equilibrium_derivative(of='rho', wrt='T', const='H')
        print("two phase")

    dpdh = drho_dT/res.dH_dT_P() # using drho/dT for constante enthalpy, is this correct?
    print("dpdh = ", dpdh)
    print("drho_dT = ", drho_dT)

    return dpdh

# 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, minL, mS, cpS, u, area):
    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] = Tin
    # y[2] = minL

    i=0

    global equil_state

    ## Using thermo
    T, equil_state = temp_thermo(y[i],P,last_equil_state=equil_state) #initial T
    m, equil_state = mass_thermo(y[i],P,last_equil_state=equil_state) #initial m

    #initial values
    dydt[i] = (minL * (hIn - y[i]) + (u * area * (y[i + 1] - T))) / m # dH/dt (eq. 2)
    dydt[i + 1] = -(u * area * (y[i + 1] - T)) / (mS * cpS) # dTs/dt from eq.3
    dmdt = dydt[i] * dpdh(y[i], P) * volumeV / n # dm/dt (holdup variation) eq. 4b
    dydt[i + 2] = minL - dmdt # mass flow out (eq.4a)


    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
        print("i = ", i)
        ## Using thermo
        T, equil_state = temp_thermo(y[i],P,last_equil_state=equil_state)
        m, equil_state = mass_thermo(y[i],P,last_equil_state=equil_state)

        # [h, TS, m]
        dydt[i] = (dydt[i-1] * (y[i - 3] - y[i]) + (u * area * (y[i + 1] - T))) /  m # dH/dt (eq.2), dydt[i-1] is the mass of the previous tank
        dydt[i + 1] = -((u * area * (y[i + 1] - T))) / (mS * cpS) # dTs/dt eq. (3)
        dmdt = dydt[i] * dpdh(y[i], P) * volumeV / n # Equation 4b
        dydt[i + 2] = dydt[i-1] - dmdt # Equation 4a

    return dydt

Ts0 = 80  # solid starting temperature (K)
y0 = np.copy(np.asarray([H0, Ts0, minL] * n)) #setting up initial condition array (H0,Ts0 and m0 at all points
y0[0] = H0 # initial fluid bed enthalpy
y0[1] = Ts0 # initial solid bed temperature
y0[2] = minL #initial mass flow        #####MINL

parameters = (n, hIn, minL, mS, cpS, u, area)

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