'''
    Suit for the Geldart A&B particles
    Calculate the U_mf U_mb rho_e mu_e epsilon_mf
    from the paper:
        "Gao X, Wang L J, Wu C, et al. Novel bubble–emulsion hydrodynamic model for 
        gas–solid bubbling fluidized beds[J]. Industrial & Engineering Chemistry Research, 
        2013, 52(31): 10835-10844."
        equation 1~4
    and the book of EMMS 
'''

from math import pow,exp

disp = 'yes'   # whether to display, yes or no

# ---------------- particle properties ---------------- #
rho_p = 4332.6# kg/m3
d_mean = 48.98e-6# sauter diameter, m, mm e-3, um e-6
F45 = 0.25# the mass fraction of particles less than 45um

# ---------------- fluid properties ---------------- #
rho_f = 0.1828# kg/m3
mu_f = 3.98e-5# N·s/m2

# ---------------- other constants ---------------- #
g = 9.81


# ---------------- start solve ---------------- #
# minimum fluidized velocity
U_mf = d_mean * d_mean * (rho_p - rho_f) * g / 1650 / mu_f
Re_p = d_mean * rho_f * U_mf / mu_f
if Re_p > 1000:
    U_mf = math.sqrt(d_mean * (rho_p - rho_f) * g / 24.5 / rho_f)
# minimum bubble velocity
U_mb = U_mf * 2300 * pow(rho_f, 0.126) * pow(mu_f , 0.523) * \
    exp(0.716 * F45) / pow(d_mean , 0.8) / pow(g , 0.954)\
        /pow(rho_p-rho_f , 0.934)

# solve two equations for :
#       epsilon_mf (minimum fluidized solid fraction) , 
#       delta_e (fluid fraction of emulsion phase)
from sympy import solve,symbols
epsilon_mf = symbols('epsilon_mf')
function1 = 1.75 * Re_p * Re_p / epsilon_mf ** 3 + 150 * (1 - epsilon_mf)\
    / epsilon_mf ** 3 * Re_p - d_mean ** 3 * rho_f * (rho_p - rho_f) * g / mu_f ** 2
solves1 = solve(function1, epsilon_mf)
# for i in solves1:
#     if i>0 && i<1:
#         epsilon_mf = i
#         break
    
# # default epsilon_mf = 0.5
epsilon_mf = float(solves1[0])

delta_e = symbols('delta_e')
constant1 = epsilon_mf ** 3 / (1 - epsilon_mf) * pow(U_mb/U_mf , 0.7)
function2 = constant1 * (1 - delta_e) - delta_e ** 3
solves2 = solve(function2 , delta_e)
# for i in solves2:
#     if i>0 && i<1:
#         delta_e = i
#         break
# # default delta_e = 0.5
delta_e = float(solves2[0])
    
alpha_p = 1 - delta_e
rho_e = rho_p * alpha_p + rho_f * delta_e
mu_e = mu_f * (1 + 2.5 * alpha_p + 10.05 * pow(alpha_p , 2) + 0.00273 * exp (16.6 * alpha_p))

if disp == 'yes':
    print("U_mb\t",'%.4f' % U_mb,'\n')
    print("U_mf\t",'%.4f' % U_mf,'\n')
    print("Re_p\t",'%.4f' % Re_p,'\n')
    print("epsilon_mf\t",'%.4f' % epsilon_mf,'\n')
    print("delta_e\t",'%.4f' % delta_e,'\n')
    print("rho_e\t",'%.4f' % rho_e,'\n')
    print("rho_p\t",'%.4f' % rho_p,'\n')
    print("rho_f\t",'%.4f' % rho_f,'\n')
    print("mu_f\t",f"{mu_f:.3E}",'\n')
    print("mu_e\t",'%.4f' % mu_e,'\n')