Steady state solution of population balance using GEKKO

user_18297833 0
Viewed 98

I'm looking to implement a model of the solution of the steady state population balance of an MSMPR crystalliser, based on the equations given here: https://pubs.acs.org/doi/abs/10.1021/acs.cgd.8b00446

For a range of values of L (say 0.5e-6 to 1000e-6), I want to solve the above for n. I've tried to do this with the GEKKO library various ways and either end up with a no solution or incorrect values. I think it's the implementation of the integral I'm getting wrong, as I can solve it in excel when M_T is known and get the expected result.

This is my current version for a value of L = 25e-6 for which n should equal 4.67e16 :

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt

m = GEKKO()

R = m.Const(value = 8.3145)
k_b0 = m.Const(value = 4.79e22) #m-3s-1
E_b = m.Const(value = 73000) #J/mol
j = m.Const(value = 0.56)
b0 = m.Const(value = 2.60)
k_g0 = m.Const(value = 13.1) #m/s
E_g = m.Const(value = 52500) #J/mol
g0 = m.Const(value = 0.87)
kv =  m.Const(value = 3.74)
rho = m.Const(value = 1280) #kg/m3 or g/L
C_0 = m.Const(value = 129.5) #g/L
Temp = m.Const(value = 25) #degC
tau = m.Const(value = 2400) #s

C_sat = m.Var()
n = m.Var()
B = m.Var()
G = m.Var()
M_T = m.Var()
L = m.Var(25e-6)
C_ml = m.Var()

m.Equation(C_sat == 15.282*m.exp(0.041*Temp))
m.Equation(n == (B/G)*m.exp(-L/(G*tau)))
m.Equation(B == k_b0*m.exp(-E_b/(R*(Temp+273.15)))*(M_T**j)*((C_0-M_T-C_sat)/C_sat)**b0)
m.Equation(G == k_g0*m.exp(-E_g/(R*(Temp+273.15)))*((C_0-M_T-C_sat)/C_sat)**g0)
m.Equation(M_T == kv*rho*m.integral(L**3*n))

m.options.SOLVER = 1
m.solve(disp=True)

Any suggestions would be appreciated.

gekko python
1 Answers

The following code may help:

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt

m = GEKKO()

R = m.Const(value = 8.3145)
k_b0 = m.Const(value = 4.79e22) #m-3s-1
E_b = m.Const(value = 73000) #J/mol
j = m.Const(value = 0.56)
b0 = m.Const(value = 2.60)
k_g0 = m.Const(value = 13.1) #m/s
E_g = m.Const(value = 52500) #J/mol
g0 = m.Const(value = 0.87)
kv =  m.Const(value = 3.74)
rho = m.Const(value = 1280) #kg/m3 or g/L
C_0 = m.Const(value = 129.5) #g/L
Temp = m.Const(value = 25) #degC
tau = m.Const(value = 2400) #s

C_sat = m.Var(40)
n = m.Var(4.67e16)
B = m.Var(1e-8)
G = m.Var(1e-9)
L = m.Const(25e-6)
M_T = m.Var(30, lb=0, ub=100)
C_ml = m.Var(20, lb=0)

m.Equation(M_T == kv*rho*6*(G*tau)**4 * (B/G))
m.Equation(C_sat == 15.282*m.exp(0.041*Temp))
m.Equation(n == (B/G)*m.exp(-L/(G*tau)))
m.Equation(B == k_b0*m.exp(-E_b/(R*(Temp+273.15)))*(M_T**j)*((C_0-M_T-C_sat)/C_sat)**b0)
m.Equation(G == k_g0*m.exp(-E_g/(R*(Temp+273.15)))*((C_0-M_T-C_sat)/C_sat)**g0)
m.Equation(C_ml == C_0 - M_T)

m.options.SOLVER = 1 # or 3
m.solve(disp=True)
user_172037030
Related