### Limestone Slurry Pipeline Optimization Problem

Design a pipeline for transporting crushed limestone from a quarry to a terminal located some distance away, using water as a transporting medium.

The limestone is crushed at the quarry, mixed with water to form a slurry, and pumped through the pipe. We would like to minimize the operating cost, which is primarily determined by the grinding power and the pumping power.

This problem is taken from James Siddall, Optimal Engineering Design, Dekker.

This assignment can be completed in groups of two. Additional guidelines on individual, collaborative, and group assignments are provided under the Expectations link.

Solution Help

from gekko import GEKKO
m = GEKKO()

# Pipeline slurry

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

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

# 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,'Rw')

# Friction factor for water (fit empirical correlation)
y = m.Intermediate(m.log(Rw),'y')
fw = m.Intermediate(-1.5919e-4*y**3+5.5535e-3*y**2-6.8029e-2*y+0.31078,'fw')

# Dimensionless numbers to obtain drag coefficient
CdRp2 = m.Intermediate(4*g*rho_w*d**3*(gamma_L-rho_w)/(3*mu**2),'CdRp2')
x = m.Intermediate(m.log(CdRp2),'x')
Cd = m.Intermediate(m.exp(0.03420*x**2-0.98327*x+6.17176),'Cd')

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

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

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

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

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

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

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

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

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

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

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

# 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.solve()
print('Optimal Power: ' + str(Pt[0]))