Мне успешно удалось визуализировать мембранный реактор постоянного тока с помощью решателяsolve_ivp. Другими словами, трубка за трубкой, обмен веществами и температурой происходит через интерфейсы. Теперь мне было интересно, как бы это выглядело, если бы я смоделировал все это в противотоке. Мои проблемы заключаются в формулировке BC и в формулировке первоначального предположения.
Исправление:
Теперь у меня есть граничные условия, а также начальное предположение, но python теперь жалуется на следующую ошибку:
Код: Выделить всё
in xsolve
reactor[j].gas.y[i] = np.array(reactor[j].Fmol[i]/reactor[j].Fmol_ges)
~~~~~~~~~~~~~~~~^^^
ValueError: setting an array element with a sequence.
Вот мой код:
Код: Выделить всё
from sub import substance, mixture
import numpy as np
from geo import reactor_volume, catalyst, membrane
from kinetics import kinetics
from math import pi
import matplotlib.pyplot as plt
import termplotlib as tpl
from scipy.integrate import solve_bvp as solve
from controller import controller
### Controller Settings
dpio = True
adiabatio = False
thermtubeio = True
thermsweepio = True
thermmemio = True
drainio = True
reacio = False
### Controller init
ctl = controller(dpio,adiabatio,thermtubeio,thermsweepio,thermmemio,drainio,reacio)
ctl.setcontroll()
## fixed variables
### Tube
Tref = 273 # reference temperature [K]
T0_t = 493 # temperature of tube inlet [K]
p0_t = 50*1e5 # pressure of tube inlet [Pa]
tube_diam = 0.016 # diameter of tube [m]
mdot_feed_t = 2.8*1e-5 # mass flow of tube inlet [kg/s]
reactor_length = 0.15 # length of reactor [m]
### Sweep
Rp = 0.1 # factor for controlling pressure in sweep
Rf = 0.1 # factor for controlling mass flow in sweep
T0_s = 400 # temperature of sweep inlet [K]
p0_s = Rp * p0_t # pressure of sweep [Pa]
outer_diam = 0.008 # outer diam only sweep
sweep_diam = tube_diam + outer_diam # diameter of sweep + tube [m]
mdot_feed_s = Rf * mdot_feed_t # mass flow of sweep inlet [kg/s]
### Catalyst
material = "CuZnO"
kat_mass = 0.0348 # mass of catalyst [kg]
kat_diam = 0.0005 # diameter of catalyst particle [m]
kat_rho_part = 1775 # density of catalyst [kg/m3]
### Substances
Q = np.array([0, 0, 0, 0, 10**(-7), 0])
M = np.array([2e-3, 28e-3, 44e-3, 32.04e-3, 18e-3, 14e-3])
y0_t = np.array([0.82, 0.04, 0.03, 0, 0, 0.11])
y0_s = np.array([1, 0.0, 0.0, 0, 0, 0.0])
Pc = [20.5, 35*1e5, 73.83*1e5, 80.9464*1e5, 220.6*1e5, 34*1e5]
Tc = [43.6, 132.92, 304.21, 512.6, 647.096, 126.2]
Tb = [79, 81.7, 194.7, 337.668, 373.2, 77.355]
acent=[-0.216, 0.048, 0.2236, 0.556, 0.038, 0.038]
name = ["H2", "CO", "CO2", "MeOH", "H2O", "N2"]
pol = [False,True,False,True,True,False]
lambda_wall = 30
### object init
sub_t, sub_s = [], []
rs = ["t","s"]
for j in range(len(rs)):
for i in range(len(M)):
exec(f"substances{j}{i}=substance(name[{i}],{i},pol[{i}],M[{i}],y0_{rs[j]}[{i}],Q[{i}],Pc=Pc[{i}],Tc=Tc[{i}],acent=acent[{i}],Tb=Tb[{i}])")
exec(f"sub_{rs[j]}.append(substances{j}{i})")
mix1,mix2 = mixture("tube", sub_t), mixture("sweep", sub_s)
gases = [mix1, mix2]
catalyst = catalyst(material,kat_diam,kat_mass,kat_rho_part)
kin = kinetics("Methanol", "Fromment")
tube, sweep = reactor_volume("tube",tube_diam,T0_t,p0_t,mdot_feed_t, mix1,catalyst,kin), reactor_volume("sweep",sweep_diam,T0_s,p0_s, mdot_feed_s,mix2)
reactor = [tube,sweep]
membrane = membrane(0,Q)
### Starting conditions
[exec(f"gases[{j}].setMolFrac(y0_{rs[j]})") for j in range(len(rs))]
[exec(f"gases[{j}].calcRho(T0_{rs[j]},p0_{rs[j]})") for j in range(len(rs))]
[exec(f"reactor[{j}].setMdot(mdot_feed_{rs[j]})") for j in range(len(rs))]
[exec(f"reactor[{j}].calcFmol(p0_{rs[j]},T0_{rs[j]})") for j in range(len(rs))]
for j in range(len(rs)):
reactor[j].setFmol(reactor[j].Fmol_ges * reactor[j].gas.y)
reactor[0].setTemp(T0_t)
reactor[1].setTemp(T0_s)
### Solve
def xsolve(t,s):
reactor[0].setPressure(s[0])
reactor[0].setFmol(s[1:7])
reactor[0].setTemp(s[7])
reactor[1].setPressure(s[8])
reactor[1].setFmol(s[9:15])
reactor[1].setTemp(s[15])
sdot = [0 for _ in range(len(s))]
for j in range(len(reactor)):
for i in range(len(M)):
reactor[j].gas.y[i] = np.array(reactor[j].Fmol[i]/reactor[j].Fmol_ges)
#reactor[j].gas.y[i]
#y_array = [[reactor[j].Fmol[i] / reactor[j].Fmol_ges for i in range(len(M))] for j in range(len(reactor))]
#y_array = [[][]]
for i in range(len(M)):
reactor[0].gas.substances[i].calcCpx(s[7])
reactor[1].gas.substances[i].calcCpx(s[15])
reactor[0].gas.substances[i].calcViscox(s[7],s[0])
reactor[1].gas.substances[i].calcViscox(s[15],s[8])
reactor[0].gas.substances[i].calcLambda(s[7])
reactor[1].gas.substances[i].calcLambda(s[15])
reactor[0].gas.calcRho(s[7],s[0])
reactor[1].gas.calcRho(s[15],s[8])
for i in range(len(reactor)):
reactor[i].gas.calcCpm()
reactor[i].gas.calcVisc()
reactor[0].gas.calcLambdam(s[7])
reactor[1].gas.calcLambdam(s[15])
reactor[0].calcVelo(s[7],s[0],inner_diam=0)
reactor[1].calcVelo(s[15],s[8],inner_diam=tube_diam)
reactor[0].kat.calcPorosity(tube_diam)
reactor[0].kat.calcRhobed()
reactor[0].calcLength()
reactor[1].length = reactor[0].length
for i in range(len(reactor)):
reactor[i].calcRe(kat_diam=kat_diam,eps=reactor[0].kat.porosity,inner_diam=tube_diam)
reactor[i].calcPr()
reactor[i].calcNu(inner_diam=tube_diam)
reactor[i].calcAlpha_wall()
reactor[0].kin.calcKeq(s[7])
reactor[0].kin.calcHr(s[7],[reactor[0].gas.substances[i].cpx for i in range(len(M))])
reactor[0].kin.calcReacrate(s[0],reactor[0].gas.y)
membrane.setInterface(reactor[1].gas.y,reactor[0].gas.y,s[8],s[0],reactor[1].alpha_wall,reactor[0].alpha_wall)
membrane.calcJ()
membrane.calcU(sweep_diam,tube_diam,lambda_wall)
sdot[0] = ctl.dpio * (-150 * reactor[0].gas.visc * (1 - reactor[0].kat.porosity)**2 * reactor[0].v / reactor[0].kat.porosity**3 / kat_diam**2 - 1.75 * (1 - reactor[0].kat.porosity) * reactor[0].gas.rho_T * reactor[0].v**2 / reactor[0].kat.porosity**3 / kat_diam)
sdot[8] = 0
for i in range(len(M)):
sdot[i+1] = ctl.reacio * reactor[0].kat.rho_bed * np.pi * tube_diam**2 * reactor[0].kin.ri[i] - membrane.J[i] * 2 * np.pi * tube_diam
sdot[i+9] = membrane.J[i] * 2 * np.pi * tube_diam
sumJicpi = ctl.drainio * sum(membrane.J[i] * reactor[0].gas.substances[i].cpx for i in range(len(membrane.J)))
sumFmcpitube = sum(reactor[0].Fmol[i] * reactor[0].gas.substances[i].cpx for i in range(len(reactor[0].Fmol)))
sumFmcpisweep = sum(reactor[1].Fmol[i] * reactor[1].gas.substances[i].cpx for i in range(len(reactor[1].Fmol)))
sdot[7] = ctl.thermtubeio * 2 * np.pi * tube_diam * (ctl.thermmemio * -membrane.memU * (reactor[0].T - reactor[1].T) - ctl.reacio* 0.5 * tube_diam* reactor[0].kat.rho_bed * (reactor[0].kin.dHr[0] * reactor[0].kin.reac_rate[0] + reactor[0].kin.dHr[1] * reactor[0].kin.reac_rate[1] + reactor[0].kin.dHr[2] * reactor[0].kin.reac_rate[2]) - sumJicpi * (reactor[0].T - Tref) ) / sumFmcpitube
sdot[15] = ctl.thermsweepio * 2 * np.pi * tube_diam * (ctl.thermmemio * membrane.memU * (reactor[0].T - reactor[1].T) + sumJicpi * (reactor[0].T - Tref) ) / sumFmcpisweep
return np.vstack(sdot)
def bc(pt,FmoltH,FmoltCO,FmoltCO2,FmoltMeOH,FmoltH2O,FmoltN,Tt,ps,FmolsH,FmolsCO,FmolsCO2,FmolsMeOH,FmolsH2O,FmolsN,Ts):
return (np.array(
reactor[0].p_ges - reactor[0].p_ges,
reactor[0].Fmol[0] - reactor[0].Fmol[0],
reactor[0].Fmol[1] - reactor[0].Fmol[1],
reactor[0].Fmol[2] - reactor[0].Fmol[2],
reactor[0].Fmol[3] - reactor[0].Fmol[3],
reactor[0].Fmol[4] - reactor[0].Fmol[4],
reactor[0].Fmol[5] - reactor[0].Fmol[5],
reactor[0].T - reactor[0].T,
reactor[1].p_ges - reactor[0].p_ges,
reactor[1].Fmol[0] - reactor[1].Fmol[0],
reactor[1].Fmol[1] - reactor[1].Fmol[1],
reactor[1].Fmol[2] - reactor[1].Fmol[2],
reactor[1].Fmol[3] - reactor[1].Fmol[3],
reactor[1].Fmol[4] - reactor[1].Fmol[4],
reactor[1].Fmol[5] - reactor[1].Fmol[5],
reactor[1].T - reactor[1].T))
t_span = np.array([0, reactor_length])
t = np.linspace(t_span[0],t_span[1],25)
s0 = np.zeros((16, t.size))
s0 = np.ones((16,t.shape[0]))
sol = solve(xsolve,bc,t,s0,tol=1e-9,verbose=2)
С уважением
Понимайте и правильно используйте решательsolve_bvp
Подробнее здесь: https://stackoverflow.com/questions/786 ... bvp-solver
