Hydrogenation Reactor

main.go

package main

import (
	"flag"
	"fmt"
	"math"

	"github.com/cpmech/gosl/la"
	"github.com/cpmech/gosl/ode"
	"github.com/cpmech/gosl/plt"
)

func main() {
	flowH2O := flag.Float64("H2O", 0, "initial flow of steam (mol/s)")
	flowH2S := flag.Float64("H2S", 0.000736, "initial flow of hydrogen sulfide (mol/s)")
	flowH2 := flag.Float64("H2", 2.455, "initial flow of hydrogen (mol/s)")
	flowT := flag.Float64("Thiophene", 0.001105, "initial flow of thiophene (mol/s)")

	flowCH4 := flag.Float64("CH4", 109.73, "initial flow of methane (mol/s)")
	flowN2 := flag.Float64("N2", 0.7365, "initial flow of nitrogen (mol/s)")
	flowC2H6 := flag.Float64("C2H6", 9.8196, "initial flow of ethane (mol/s)")
	flowCO2 := flag.Float64("CO2", 2.4549, "initial flow of CO2 (mol/s)")

	ρ := flag.Float64("density", 11.35, "gas density (kg/m^3)")
	T := flag.Float64("T", 673.15, "initial reactor temperature (K)")
	P := flag.Float64("P", 2450, "initial reactor pressure (kPa)")

	ϕ := flag.Float64("voidage", 0.45, "bed voidage (ϕ)")
	μ := flag.Float64("viscosity", 0.0000215, "gas viscosity (μ)")
	Dp := flag.Float64("Dp", 0.0025, "particle diameter (m)")
	ρb := flag.Float64("catalyst-density", 590, "catalyst-density (kg/m^3)")

	plot := flag.Bool("plot", false, "stops plotting of graphs")

	D := flag.Float64("d", 2.0197, "bed diameter (m)")

	flag.Parse()
	ρc := *ρb / (1.0 - *ϕ)
	area := math.Pow(*D, 2) * math.Pi / 4
	F0 := *flowH2 + *flowT + *flowH2S + *flowH2O + *flowC2H6 + *flowCH4 + *flowCO2 + *flowN2
	W := 1500.0

	ODEs := func(f la.Vector, h, x float64, y la.Vector) {
		totalFlow := y[0] + y[1] + y[2] + y[3] + *flowC2H6 + *flowCH4 + *flowCO2 + *flowN2
		flows := map[string]float64{
			"H2O":  y[0],
			"H2":   y[1],
			"T":    y[2],
			"H2S":  y[3],
			"C2H6": *flowC2H6,
			"CH4":  *flowCH4,
			"CO2":  *flowCO2,
			"N2":   *flowN2,
		}
		partials := map[string]float64{}
		for compound, flow := range flows {
			partials[compound] = flow / totalFlow * y[5]
		}
		G := (48.11*flows["T"] + 34.1*flows["H2S"] + 2.016*flows["H2"] + 18.015*flows["H2O"] + 16.04*flows["CH4"] + 28.013*flows["N2"] + 30.07*flows["C2H6"] + 44.01*flows["CO2"]) / 1000 / area
		beta := β(*ϕ, G, *Dp, *μ, *ρ)
		alpha := α(beta, area, ρc, *ϕ, *P)

		f[0] = dFH2OdW(y[4], partials)
		f[1] = dFH2dW(y[4], partials)
		f[2] = dFTdW(y[4], partials)
		f[3] = dFH2SdW(y[4], partials)
		f[4] = dTdW(y[4], flows, partials)
		f[5] = dPdW(alpha, y[5], *P, y[4], *T, totalFlow, F0)
	}

	config := ode.NewConfig("radau5", "", nil)
	config.SetStepOut(true, nil)

	parameters := []float64{*flowH2O,
		*flowH2,
		*flowT,
		*flowH2S,
		*T, *P}

	solver := ode.NewSolver(len(parameters), config, ODEs, nil, nil)
	defer solver.Free()
	solver.Solve(la.NewVectorSlice(parameters), 0, W)

	wValues := solver.Out.GetStepX()
	yValues := solver.Out.GetStepYtableT()

	conversion := 1 - yValues[2][len(wValues)-1]/yValues[2][0]
	pressureDrop := *P - yValues[5][len(wValues)-1]

	fmt.Printf("conversion: %.2f; pressure drop (kPa): %.4f; outlet temperature %2f (K)\n", conversion, pressureDrop, yValues[4][len(wValues)-1])
	flows := []float64{parameters[0],
		parameters[1],
		parameters[2],
		parameters[3],
	}

	fmt.Printf("flows (mol/s); H2O: %.2f; H2: %.2f; T: %.2f; H2S: %.2f; total catalyst: %.2f kg\n", flows[0], flows[1], flows[2], flows[3], W)
	var conversions []float64

	for _, v := range yValues[2] {
		conversions = append(conversions, 1-v/yValues[2][0])
	}

	if !*plot {
		return
	}

	plt.Subplot(1, 3, 1)
	plt.Plot(wValues, conversions, nil)
	plt.Grid(nil)
	plt.SetLabels("Catalyst (kg)", "Thiophene Conversion", nil)

	plt.Subplot(1, 3, 2)
	plt.Plot(wValues, yValues[5], nil)
	plt.SetTicksNormal()
	plt.Grid(nil)
	plt.SetLabels("Catalyst (kg)", "Presssure (kPa)", nil)

	plt.Subplot(1, 3, 3)
	plt.Plot(wValues, yValues[4], nil)
	plt.SetTicksNormal()
	plt.Grid(nil)
	plt.SetLabels("Catalyst (kg)", "T (K)", nil)

	plt.Show()
}

ode.go

package main

import (
	"math"

	"github.com/ewancook/reactor/thermo"
)

func dFH2OdW(T float64, partials map[string]float64) float64 {
	return reaction(T, partials)
}

func dFTdW(T float64, partials map[string]float64) float64 {
	return -reaction(T, partials)
}

func dFH2SdW(T float64, partials map[string]float64) float64 {
	return reaction(T, partials)
}

func dFH2dW(T float64, partials map[string]float64) float64 {
	return -reaction(T, partials)
}

func dTdW(T float64, flows, partials map[string]float64) float64 {
	var denominator float64
	for compound, flow := range flows {
		denominator += thermo.SpecificHeat(compound, T) * flow
	}
	heats := reaction(T, partials) * reactionEnthalpy(T)
	return -heats * 1000 / denominator
}

func dPdW(alpha, P, P0, T, T0, F, F0 float64) float64 {
	return -alpha / 2 * P0 / (P / P0) * (T / T0) * (F / F0)
}

func α(beta, area, ρc, ϕ, P0 float64) float64 {
	return 2 * beta / (area * ρc * (1 - ϕ) * P0 * 1000)
}

func β(ϕ, G, Dp, μ, ρg float64) float64 {
	return (G * (1 - ϕ) / (ρg * Dp * math.Pow(ϕ, 3))) * (1.75*G + 150*(1-ϕ)*μ/Dp)
}

reactions.go

package main

import (
	"math"

	. "github.com/ewancook/reactor/thermo"
)

const R = 8.314

func reactionEnthalpy(T float64) float64 {
	return Enthalpy("H2O", T) + Enthalpy("H2S", T) - Enthalpy("H2", T) - 120
}

func k(T float64) float64 {
	return 1500 * math.Exp(-52000/R/T)
}

func KTH(T float64) float64 {
	return 8.42 / 100 * math.Exp(189/R/T)
}

func KH2(T float64) float64 {
	return 52210 / 100 * math.Exp(930/R/T)
}

func KH2S(T float64) float64 {
	return 19.64 / 100 * math.Exp(878/R/T)
}

func reaction(T float64, partials map[string]float64) float64 {
	return k(T) / 60 * KTH(T) * KH2(T) * partials["T"] * partials["H2"] / math.Pow(1+math.Sqrt(KH2(T))*math.Sqrt(partials["H2"]), 2) / (1 + KTH(T)*partials["T"] + KH2S(T)*partials["H2S"]/KH2(T)/partials["H2"])
}