Hi there smart people!
I am trying to implement a 1D, steady-state, real gas (compressibility factor) pipe flow model in Python using Pyomo + SCIP. It basically amounts to solving a DAE system. The formulation is an adopted version from chapter 10 in Koch, T.; Hiller, B.; Pfetsch, M.E.; Schewe, L. Evaluating Gas Network Capacities; Series on Optimization, MOS-SIAM, 2015.
However, I am encountering several problems:
- The problem seems to be numerically sensitive to a combination of the discretization step length and input parameters (mass flow, pipe length).
- Does not solve for any other model but ideal gas.
With a suitable discretization, and an ideal gas law, I get a result that seems reasonable (see example). In all other cases it turns out to be infeasible. I may be overlooking something here, but I do not see it. Therefore, if anyone is inclined to try and help me out here, I would be thankful.
The example below should produce a valid output.
Edit: I realized I had one false constraint in there belonging to another model. The energy equation works now. However, the problems mentioned above remain.
from pyomo.dae import *
from pyomo.environ import *
import matplotlib.pyplot as plt
from math import pi
from dataclasses import dataclass
@dataclass
class ShomateParameters:
A: float
B: float
C: float
D: float
E: float
def specific_isobaric_heat_capacity(self, temperature):
# in J/(mol*K)
return self.A + self.B * (temperature/1000) + self.C * (temperature/1000)**2 + self.D * (temperature/1000)**3 + self.E/(temperature/1000)**2
def plot(self, temperature_start, temperature_end, points_to_mark=None):
assert temperature_start <= temperature_end, "temperature_start <= temperature_end must hold."
temperatures = [i for i in range(temperature_start, temperature_end+1)]
values = [self.specific_isobaric_heat_capacity(temp) for temp in temperatures]
fig, ax = plt.subplots()
ax.plot(temperatures, values)
if points_to_mark is not None:
ax.scatter(points_to_mark, [self.specific_isobaric_heat_capacity(temp) for temp in points_to_mark])
ax.set(xlabel='temperature [K]', ylabel='specific isobaric heat capacity [J/(mol*K)]',
title='Shomate equation:\n A + B*T + C*T^2 + D*T^3 + E/T^2')
ax.grid()
plt.show()
@dataclass
class Species:
MOLAR_MASS: float # kg/mol
CRITICAL_TEMPERATURE: float # Kelvin
CRITICAL_PRESSURE: float # Pa
DYNAMIC_VISCOSITY: float # Pa*s
SHOMATE_PARAMETERS: ShomateParameters
METHANE = Species(MOLAR_MASS=0.016043,
CRITICAL_TEMPERATURE=190.56,
CRITICAL_PRESSURE=4599000,
DYNAMIC_VISCOSITY=1.0245e-5,
SHOMATE_PARAMETERS=ShomateParameters(
A=-0.703029,
B=108.4773,
C=-42.52157,
D=5.862788,
E=0.678565))
# select gas species
gas = METHANE
# select equation of state ('ideal', 'AGA' or 'Papay')
formula = 'ideal'
PIPE_LENGTH = 24000 # meter
start = 0 # meter
end = start + PIPE_LENGTH
MASS_FLOW = 350 # kg/s
PIPE_SLOPE = 0.0
PIPE_DIAMETER = 1.0 # meter
PIPE_INNER_ROUGHNESS = 6e-5 # 15e-6 # meter 6e-6 # meter
# gravitational acceleration
G = 9.81 # meter**2/s**2
# gas temperature at entry
TEMPERATURE = 283.15
# temperature ambient soil
TEMPERATURE_SOIL = 283.15 # Kelvin
# gas pressure at entry
PRESSURE = 4.2e6 # Pa
GAS_CONSTANT = 8.314 # J/(mol*K)
print(gas.SHOMATE_PARAMETERS)
print(gas.SHOMATE_PARAMETERS.specific_isobaric_heat_capacity(TEMPERATURE))
gas.SHOMATE_PARAMETERS.plot(273, 400, points_to_mark=[TEMPERATURE])
##################################################################################
# Variable bounds
##################################################################################
pressure_bounds = (0, 1e7) # Pa
temperature_bounds = (0, 650) # Kelvin
density_bounds = (0, 100)
idealMolarIsobaricHeatCapacityBounds = (0, 200)
correctionIdealMolarIsobaricHeatCapacityBounds = (-250, 250)
velocity_bounds = (0, 300)
##################################################################################
# Create model
##################################################################################
m = ConcreteModel()
##################################################################################
# Parameters
##################################################################################
m.criticalPressure = Param(initialize=gas.CRITICAL_PRESSURE)
m.criticalTemperature = Param(initialize=gas.CRITICAL_TEMPERATURE)
m.molarMass = Param(initialize=gas.MOLAR_MASS)
m.dynamicViscosity = Param(initialize=gas.DYNAMIC_VISCOSITY)
m.gravitationalAcceleration = Param(initialize=G)
m.pipeSlope = Param(initialize=PIPE_SLOPE)
m.innerPipeRoughness = Param(initialize=PIPE_INNER_ROUGHNESS)
m.c_HT = Param(initialize=2)
m.pi = Param(initialize=pi)
m.temperatureSoil = Param(initialize=TEMPERATURE_SOIL)
m.gasConstantR = Param(initialize=GAS_CONSTANT)
m.massFlow = Param(initialize=MASS_FLOW)
m.pipeDiameter = Param(initialize=PIPE_DIAMETER)
m.crossSectionalArea = Param(initialize=m.pi * m.pipeDiameter**2 / 4)
m.alpha = Param(initialize=3.52)
m.beta = Param(initialize=-2.26)
m.gamma = Param(initialize=0.274)
m.delta = Param(initialize=-1.878)
m.e = Param(initialize=2.2)
m.d = Param(initialize=2.2)
##################################################################################
# Variables
##################################################################################
m.x = ContinuousSet(bounds=(start, end))
m.pressure = Var(m.x, bounds=pressure_bounds) #
m.dpressuredx = DerivativeVar(m.pressure, wrt=m.x, initialize=0, bounds=(-100, 100))
m.temperature = Var(m.x, bounds=temperature_bounds) #
m.dtemperaturedx = DerivativeVar(m.temperature, wrt=m.x, initialize=0, bounds=(-100, 100))
m.density = Var(m.x, bounds=density_bounds)
m.ddensitydx = DerivativeVar(m.density, wrt=m.x, initialize=0, bounds=(-100, 100))
m.z = Var(m.x, bounds=(-10, 10))
m.specificIsobaricHeatCapacity = Var(m.x)
m.idealMolarIsobaricHeatCapacity = Var(m.x, bounds=idealMolarIsobaricHeatCapacityBounds)
m.correctionIdealMolarIsobaricHeatCapacity = Var(m.x, bounds=correctionIdealMolarIsobaricHeatCapacityBounds)
m.mue_jt = Var(bounds=(-100, 100))
m.velocity = Var(m.x, bounds=velocity_bounds)
m.phiVar = Var()
##################################################################################
# Constraint rules
##################################################################################
# compressibility factor z and its derivatives; (pV/(nRT)=z
def z(p,
T,
p_c,
T_c,
formula=None):
p_r = p/p_c
T_r = T/T_c
if formula is None:
raise ValueError("formula must be equal to 'AGA' or 'Papay' or 'ideal'")
elif formula == 'AGA':
return 1 + 0.257 * p_r - 0.533 * p_r/T_r
elif formula == 'Papay':
return 1-3.52 * p_r * exp(-2.26 * T_r) + 0.247 * p_r**2 * exp(-1.878 * T_r)
elif formula == 'ideal':
return 1
else:
raise ValueError("formula must be equal to 'AGA' or 'Papay' or 'ideal'")
def dzdT(p,
T,
p_c,
T_c,
formula=None):
p_r = p/p_c
T_r = T/T_c
if formula is None:
raise ValueError("formula must be equal to 'AGA' or 'Papay'")
elif formula == 'AGA':
return 0.533 * p/p_c * T_c * 1/T**2
elif formula == 'Papay':
return 3.52 * p_r * (2.26/T_c) * exp(-2.26 * T_r) + 0.247 * p_r**2 * (-1.878/T_c) * exp(-1.878 * T_r)
elif formula == 'ideal':
return 0
else:
raise ValueError("formula must be equal to 'AGA' or 'Papay' or 'ideal'")
def dzdp(p,
T,
p_c,
T_c,
formula=None):
p_r = p/p_c
T_r = T/T_c
if formula is None:
raise ValueError("formula must be equal to 'AGA' or 'Papay' or 'ideal'")
elif formula == 'AGA':
return 0.257 * 1/p_c - 0.533 * (1/p_c)/T_r
elif formula == 'Papay':
return -3.52 * 1/p_c * exp(-2.26 * T_r) + 0.274 * 1/(p_c**2) * 2 * p * exp(-1.878 * T_r)
elif formula == 'ideal':
return 0
else:
raise ValueError("formula must be equal to 'AGA' or 'Papay' or 'ideal'")
def make_c_compr(formula):
assert formula == 'AGA' or formula == 'Papay' or formula == 'ideal'
def _c_compr(z_var,
p,
T,
p_c,
T_c):
return z_var - z(p, T, p_c, T_c, formula=formula)
return _c_compr
_c_compr = make_c_compr(formula)
def _c_compr_rule(m, x):
return 0 == _c_compr(m.z[x],
m.pressure[x],
m.temperature[x],
m.criticalPressure,
m.criticalTemperature)
m.c_compr = Constraint(m.x, rule=_c_compr_rule)
# specific isobaric heat capacity: ideal + correction term
def _c_mhc_real(molarMass,
specificIsobaricHeatCapacity,
idealMolarIsobaricHeatCapacity,
correctionIdealMolarIsobaricHeatCapacity):
return molarMass * specificIsobaricHeatCapacity - (idealMolarIsobaricHeatCapacity +
correctionIdealMolarIsobaricHeatCapacity)
def _c_mhc_real_rule(m, x):
return 0 == _c_mhc_real(m.molarMass,
m.specificIsobaricHeatCapacity[x],
m.idealMolarIsobaricHeatCapacity[x],
m.correctionIdealMolarIsobaricHeatCapacity[x])
m.c_mhc_real = Constraint(m.x, rule=_c_mhc_real_rule)
# _c_mhc_ideal_Shomate
def _c_mhc_ideal_Shomate(idealMolarIsobaricHeatCapacity, A, B, C, D, E, T):
return idealMolarIsobaricHeatCapacity - (A + B * (T/1000) + C * (T/1000)**2 + D * (T/1000)**3 + E/(T/1000)**2)
def _c_mhc_ideal_Shomate_rule(m, x):
return 0 == _c_mhc_ideal_Shomate(m.idealMolarIsobaricHeatCapacity[x],
gas.SHOMATE_PARAMETERS.A,
gas.SHOMATE_PARAMETERS.B,
gas.SHOMATE_PARAMETERS.C,
gas.SHOMATE_PARAMETERS.D,
gas.SHOMATE_PARAMETERS.E,
m.temperature[x])
m.c_mhc_ideal_Shomate = Constraint(m.x, rule=_c_mhc_ideal_Shomate_rule)
# _c_mhc_corr
def make_c_mhc_corr(formula):
assert formula == 'AGA' or formula == 'Papay' or formula == 'ideal'
if formula == 'AGA':
def _c_mhc_corr(correctionIdealMolarIsobaricHeatCapacity, alpha, beta, gamma, delta, p, T, p_c, T_c, R):
return correctionIdealMolarIsobaricHeatCapacity
elif formula == 'Papay':
def _c_mhc_corr(correctionIdealMolarIsobaricHeatCapacity, alpha, beta, gamma, delta, p, T, p_c, T_c, R):
# m.alpha = 3.52
# m.beta = -2.26
# m.gamma = 0.274
# m.delta = -1.878
p_r = p/p_c
T_r = T/T_c
return correctionIdealMolarIsobaricHeatCapacity + R * (
(gamma * delta + 0.5 * gamma * delta**2 * T_r) * p_r**2 * T_r * exp(delta * T_r) -
(2 * alpha * beta + alpha * beta**2 * T_r) * p_r * T_r * exp(beta * T_r))
elif formula == 'ideal':
def _c_mhc_corr(correctionIdealMolarIsobaricHeatCapacity, alpha, beta, gamma, delta, p, T, p_c, T_c, R):
return correctionIdealMolarIsobaricHeatCapacity
return _c_mhc_corr
_c_mhc_corr = make_c_mhc_corr(formula)
def _c_mhc_corr_rule(m, x):
return 0 == _c_mhc_corr(m.correctionIdealMolarIsobaricHeatCapacity[x],
m.alpha,
m.beta,
m.gamma,
m.delta,
m.pressure[x],
m.temperature[x],
m.criticalPressure,
m.criticalTemperature,
m.gasConstantR)
m.c_mhc_corr = Constraint(m.x, rule=_c_mhc_corr_rule)
# equation of state
def _c_eos(p, T, rho, molarMass, R, z):
return rho * z * R * T - p * molarMass
def _c_eos_rule(m, x):
return 0 == _c_eos(m.pressure[x],
m.temperature[x],
m.density[x],
m.molarMass,
m.gasConstantR,
m.z[x])
m.c_eos = Constraint(m.x, rule=_c_eos_rule)
# flow velocity equation
def _c_vel_flow(q, v, rho, A):
return A * rho * v - q
def _c_vel_flow_rule(m, x):
return 0 == _c_vel_flow(m.massFlow,
m.velocity[x],
m.density[x],
m.crossSectionalArea)
m.c_vel_flow = Constraint(m.x, rule=_c_vel_flow_rule)
# a smooth reformulation of the flow term with friction: lambda(q)|q|q (=phi)
def _c_friction(phi, q, k, D, e, d, A, eta):
beta = k/(3.71 * D)
lambda_slant = 1/(2*log10(beta))**2
alpha = 2.51 * A * eta / D
delta = 2 * alpha/(beta*log(10))
b = 2 * delta
c = (log(beta) + 1) * delta**2 - (e**2 / 2)
root1 = sqrt(q**2 + e**2)
root2 = sqrt(q**2 + d**2)
return phi - lambda_slant * (root1 + b + c/root2) * q
def _c_friction_rule(m):
return 0 == _c_friction(m.phiVar,
m.massFlow,
m.innerPipeRoughness,
m.pipeDiameter,
m.e,
m.d,
m.crossSectionalArea,
m.dynamicViscosity)
m.c_friction = Constraint(rule=_c_friction_rule)
# energy balance
def _diffeq_energy(q, specificIsobaricHeatCapacity, dTdx, T, rho, z, dzdT, dpdx, g, s, pi, D, c_HT, T_soil):
return q * specificIsobaricHeatCapacity * dTdx - (q * T / (rho * z) * dzdT * dpdx) + (q * g * s) + (pi * D * c_HT * (T - T_soil))
def _diffeq_energy_rule(m, x):
# if x == start:
# return Constraint.Skip
return 0 == _diffeq_energy(m.massFlow,
m.specificIsobaricHeatCapacity[x],
m.dtemperaturedx[x],
m.temperature[x],
m.density[x],
m.z[x],
dzdT(m.pressure[x],
m.temperature[x],
m.criticalPressure,
m.criticalTemperature,
formula=formula),
m.dpressuredx[x],
m.gravitationalAcceleration,
m.pipeSlope,
m.pi,
m.pipeDiameter,
m.c_HT,
m.temperatureSoil)
m.diffeq_energy = Constraint(m.x, rule=_diffeq_energy_rule)
# momentum balance
def _diffeq_momentum(rho, dpdx, q, A, drhodx, g, s, phi, D):
return rho * dpdx - q**2 / (A**2) * drhodx / rho + g * rho**2 * s + phi / (2 * A**2 * D)
def _diffeq_momentum_rule(m, x):
# if x == start:
# return Constraint.Skip
return 0 == _diffeq_momentum(m.density[x],
m.dpressuredx[x],
m.massFlow,
m.crossSectionalArea,
m.ddensitydx[x],
m.gravitationalAcceleration,
m.pipeSlope,
m.phiVar,
m.pipeDiameter)
m.diffeq_momentum = Constraint(m.x, rule=_diffeq_momentum_rule)
##################################################################################
# Discretization
##################################################################################
discretizer = TransformationFactory('dae.finite_difference')
discretizer.apply_to(m, nfe=60, wrt=m.x, scheme='BACKWARD')
##################################################################################
# Initial conditions
##################################################################################
m.pressure[start].fix(PRESSURE)
m.temperature[start].fix(TEMPERATURE)
##################################################################################
# Objective
##################################################################################
# constant
m.obj = Objective(expr=1)
m.pprint()
##################################################################################
# Solve
##################################################################################
solver = SolverFactory('scip')
# solver = SolverFactory('scip', executable="C:/Users/t.triesch/Desktop/scipampl-7.0.0.win.x86_64.intel.opt.spx2.exe")
results = solver.solve(m, tee=True, logfile="pipe.log")
##################################################################################
# Plot
##################################################################################
distance = [i/1000 for i in list(m.x)]
p = [value(m.pressure[x])/1e6 for x in m.x]
t = [value(m.temperature[x]) for x in m.x]
rho = [value(m.density[x]) for x in m.x]
v = [value(m.velocity[x]) for x in m.x]
fig = plt.figure()
gs = fig.add_gridspec(4, hspace=0.5)
axs = gs.subplots(sharex=True)
fig.suptitle('p[start] = {0} [MPa], p[end] = {1} [MPa],\n T[start]= {2} [K],\n massFlow[:]= {3} [kg/s],\n total length: {4} m'.format(
p[0], p[-1], t[0], m.massFlow.value, PIPE_LENGTH))
axs[0].plot(distance, p, '-')
axs[0].set(ylabel='p [MPa]')
axs[0].set_ylim([0, 10])
axs[0].grid()
axs[0].set_yticks([i for i in range(0, 11)])
axs[1].plot(distance, t, '-')
axs[1].set(ylabel='T [K]')
axs[1].set_ylim([250, 350])
axs[1].grid()
axs[2].plot(distance, rho, '-')
axs[2].set(ylabel='rho [kg/m^3]')
axs[2].grid()
axs[3].plot(distance, v, '-')
axs[3].set(ylabel='v [m/s]')
axs[3].grid()
for ax in axs.flat:
ax.set(xlabel='distance [km]')
plt.show()
