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