import numpy as np
from gekko import GEKKO
m = GEKKO()

# Slurry Pipeline

# Constants
L = 15*5280 # Length of pipeline (ft)
W = 12.67   # Massflow of limestone (lbm/sec)
a = 0.01    # Average lump-size before grinding (ft)
pi = 3.1415927

rho_w = 62.428 # Density of water (lbm/ft^3)
mu = 7.392e-4  # Viscosity of water (lbm/ft/sec)
g = 32.174     # Gravitational constant (ft/sec^2)
# Conversion between lbm and lbf (lbm ft / lbf / s^2)
gc = 32.174049 
# Density of limestone (lbm/ft^3)
gamma_L = 168.5 

# Variables
# Average flow velocity (ft/sec)
V = m.Var(value=10, lb=1, ub=20, name='V')
# Volumetric concentration of slurry
#   (Vol limestone/Vol Total)
c = m.Var(value=0.2, lb=0.01, ub=0.4, name='c') 
# Internal pipe diameter (ft)
Dpipe = m.Var(value=0.4, lb=0.01, ub=0.5, name='Dpipe') 
# Average particle size after grinding (ft)
d = m.Var(value=0.008, lb=0.0008, name='d') 
Pt = m.Var(value=1, name='Pt')

# Intermediates
# Reynolds number
Rw = m.Intermediate(rho_w*V*Dpipe/mu)

# Friction factor
ffp = True
if ffp:
    y = m.Intermediate(m.log(Rw))
    fw = m.Intermediate(-1.5919e-4*y**3+5.5535e-3*y**2\
                        -6.8029e-2*y+0.31078)
else:
    fw1 = 0.3164/(Rw**0.25)
    fw2 = 0.0032+0.221*Rw**-0.237
    fw = m.if3(Rw-1e5,fw1,fw2)

# Dimensionless numbers to obtain drag coefficient
cspline = True
CdRp2 = m.Intermediate(4*g*rho_w*d**3*\
                       (gamma_L-rho_w)/(3*mu**2))
if cspline:
    import pandas as pd
    course = 'http://apmonitor.com/me575/'
    url = 'index.php/Main/LimestoneSlurry'
    data = pd.read_html(course+url)[0] # read data
    lnCdRp2,lnCd,Cd = m.Array(m.Var,3,value=2)
    m.Equation(m.exp(lnCdRp2)==CdRp2)
    m.Equation(m.exp(lnCd)==Cd)
    m.cspline(lnCdRp2,lnCd,
              np.log(data['CdRp2'].values),
              np.log(data['Cd'].values),False)
else:
    x = m.Intermediate(m.log(CdRp2))
    Cd = m.Intermediate(m.exp(0.03420*x**2\
                              -0.98327*x+6.17176))

# Slurry density
rho = m.Intermediate((1-c)*rho_w+c*gamma_L)

# Limestone specific gravity
S = m.Intermediate(gamma_L/rho_w)

# Slurry friction factor
r = m.Intermediate(rho_w/rho)
t = m.Intermediate(g*Dpipe*(S-1)/(V**2*m.sqrt(Cd)))
f = m.Intermediate(fw*(r+150*c*r*t**1.5))

# Pressure drop (lbf / ft^2)
dp = m.Intermediate(f*rho*L*V**2/(2*Dpipe*gc))

# Pipe cross-sectional area (ft^2)
Area = m.Intermediate(pi*Dpipe**2/4)

# Slurry flow rate (ft^3/sec)
Q = m.Intermediate(Area*V)

# Slurry mass flow rate (lbm/sec)
mdot = m.Intermediate(rho*Area*V)

# Limestone mass flow rate (lbm/sec)
Q_L = m.Intermediate(Q*c) # ft^3/sec
# lbm/sec = (lbm/ft^3) * (ft^3/sec)
mdot_L = m.Intermediate(gamma_L*Q_L) 

# Friction power loss (ft * lbf / sec)
Pf = m.Intermediate(dp*Q)

# Grinding power (ft*lbf/sec)
Pg = m.Intermediate(218*W*(1/m.sqrt(d)-1/m.sqrt(a)))

# Critical velocity
Vc = m.Intermediate((40*g*c*(S-1)*Dpipe/m.sqrt(Cd))**0.5)

# Equations
# Total power (grinding + pumping)
m.Equation(Pt==Pg+Pf)

# Mass flow of limestone
m.Equation(mdot_L==W)

# Velocity must be greater than 
#   the critical velocity constraint
m.Equation(V>Vc)

m.Minimize(Pt)

m.options.SOLVER=1
m.solve()
print('Optimal Power: ' + str(Pt[0]))