Scipy_integrate и Sympy для решения BVP ⇐ Python

[Anonymous]

У меня есть небольшая проблема, а именно то, что у меня 8 дифференциалов (в противоточном теплообменнике с массообменом через стенку), но я могу вызвать только каждый отдельный дифференциал. Проблема здесь в том, что в уравнениях есть переменные, которые зависят от других дифференциалов. Например, в моем массовом транспорте кинетика реакции зависит от температуры, а это другой дифференциал. Моя проблема в том, что я не знаю, как одновременно использовать следующие дифференциалы в решателе.
Я использую следующий решатель для решения задачи с двумя граничными значениями: https:// github.com/DavidPoves/BVP-Solver/blob/main/BVPInterface.py
Это мой код, и, может быть, кто-нибудь знает, как получить все дифференциалы за один вызов решателя?
from math import pi
import numpy as np
import BVPInterface as BVPInterface

Tref = 273 # reference temperature [K]
T0_t = 493 # temperature of tube inlet [K]
p0_t = 50*1e5 # pressure of tube inlet [Pa]
Rp = 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]
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])
tube_diam = 1
kat_diam = 1
reactor[0].kat.rho_bed = 1
reactor[0].kin.ri[0] = 1
reactor[0].kin.ri[1] = 2
reactor[0].kin.ri[2] = 3
reactor[0].kin.ri[3] = 4
reactor[0].kin.ri[4] = 5
reactor[0].kin.ri[5] = 6

c = 2 * np.pi * tube_diam

### Momentumtransport Equation
b = -150 / kat_diam**2
a = - 1.75 / kat_diam
dp_dL = 'b + a'

### Masstransport Equation
### H2
j1 = c
r1 = reactor[0].kat.rho_bed * np.pi * tube_diam**2 * reactor[0].kin.ri[0]
dn1H2_dL = 'r1'
dn2H2_dL = 'j1'
### CO
j2 = c
r2 = reactor[0].kat.rho_bed * np.pi * tube_diam**2 * reactor[0].kin.ri[1]
dn1CO_dL = 'r2'
dn2CO_dL = 'j2'
### CO2
j3 = c
r3 = reactor[0].kat.rho_bed * np.pi * tube_diam**2 * reactor[0].kin.ri[2]
dn1CO2_dL = 'r3'
dn2CO2_dL = 'j3'
### MeOH
j4 = c
r4 = reactor[0].kat.rho_bed * np.pi * tube_diam**2 * reactor[0].kin.ri[3]
dn1MeOH_dL= 'r4'
dn2MeOH_dL= 'j4'
### H2O
r5 = reactor[0].kat.rho_bed * np.pi * tube_diam**2 * reactor[0].kin.ri[4]
dn1H2O_dL = 'r5'
dn2H2O_dL = 'j5'
### N2
j6 = c
r6 = reactor[0].kat.rho_bed * np.pi * tube_diam**2 * reactor[0].kin.ri[5]
dn1N2_dL = 'r6'
dn2N2_dL = 'j6'

### Energytransport Equation
U = 0.5
sumJicpi = 0.5
sumFmcpitube = 0.5
sumFmcpisweep = 0.5
sumReac = 0.5

dT1_dL = 'c * ( -U * (T0_t - T0_s) - sumReac - sumJicpi * (T0_t - Tref) ) / sumFmcpitube'
dT2_dL = 'c * ( U * (T0_t - T0_s) + sumJicpi * (T0_t - Tref)) / sumFmcpisweep'

bvp_int = BVPInterface.BVPInterface()

### Define a system Momentum
syst = ['dp_dL', '0*dp_dL'] # This is the dy/dt vector
bvp_int.define_system(syst,functions=('p0_t', 'p0_s'),
vars_dict={'a': a,'b': b, 'dp_dL': dp_dL})
bvp_int.load_bcs({'a': {'1': p0_t}, 'b': {'2': p0_s}})
bvp_int.load_mesh(np.linspace(0, 0.15, 20))
sol = bvp_int.solve(tol=1e-9,verbose=2)

### Define a system H2
syst = ['dn1H2_dL-dn2H2_dL', '-dn2H2_dL'] # This is the dy/dt vector
bvp_int.define_system(syst,functions=( 'y0_t[0]', 'y0_s[0]'),
vars_dict={ 'r1':r1,'j1':j1,'dn1H2_dL': dn1H2_dL, 'dn2H2_dL':dn2H2_dL})
bvp_int.load_bcs({'a': {'1': y0_t[0]}, 'b': {'2': y0_s[0]}})
bvp_int.load_mesh(np.linspace(0, 0.15, 20))
sol = bvp_int.solve(tol=1e-9,verbose=2)

### Define a system CO
syst = ['dn1CO_dL-dn2CO_dL', '-dn2CO_dL'] # This is the dy/dt vector
bvp_int.define_system(syst,functions=( 'y0_t[1]', 'y0_s[1]'),
vars_dict={ 'r2':r2,'j2':j2,'dn1CO_dL': dn1CO_dL, 'dn2CO_dL':dn2CO_dL})
bvp_int.load_bcs({'a': {'1': y0_t[1]}, 'b': {'2': y0_s[1]}})
bvp_int.load_mesh(np.linspace(0, 0.15, 20))
sol = bvp_int.solve(tol=1e-9,verbose=2)

### Define a system CO2
syst = ['dn1CO2_dL-dn2CO2_dL', '-dn2CO2_dL'] # This is the dy/dt vector
bvp_int.define_system(syst,functions=( 'y0_t[2]', 'y0_s[2]'),
vars_dict={ 'r3':r3,'j3':j3,'dn1CO2_dL': dn1CO2_dL, 'dn2CO2_dL':dn2CO2_dL})
bvp_int.load_bcs({'a': {'1': y0_t[2]}, 'b': {'2': y0_s[2]}})
bvp_int.load_mesh(np.linspace(0, 0.15, 20))
sol = bvp_int.solve(tol=1e-9,verbose=2)

### Define a system MeOH
syst = ['dn1MeOH_dL-dn2MeOH_dL', '-dn2MeOH_dL'] # This is the dy/dt vector
bvp_int.define_system(syst,functions=( 'y0_t[3]', 'y0_s[3]'),
vars_dict={ 'r4':r4,'j4':j4,'dn1MeOH_dL': dn1MeOH_dL, 'dn2MeOH_dL':dn2MeOH_dL})
bvp_int.load_bcs({'a': {'1': y0_t[3]}, 'b': {'2': y0_s[3]}})
bvp_int.load_mesh(np.linspace(0, 0.15, 20))
sol = bvp_int.solve(tol=1e-9,verbose=2)

### Define a system H2O
syst = ['dn1H2O_dL-dn2H2O_dL', '-dn2H2O_dL'] # This is the dy/dt vector
bvp_int.define_system(syst,functions=( 'y0_t[4]', 'y0_s[4]'),
vars_dict={ 'r5':r5,'j5':j5,'dn1H2O_dL': dn1H2O_dL, 'dn2H2O_dL':dn2H2O_dL})
bvp_int.load_bcs({'a': {'1': y0_t[4]}, 'b': {'2': y0_s[4]}})
bvp_int.load_mesh(np.linspace(0, 0.15, 20))
sol = bvp_int.solve(tol=1e-9,verbose=2)

### Define a system N2
syst = ['dn1N2_dL-dn2N2_dL', '-dn2N2_dL'] # This is the dy/dt vector
bvp_int.define_system(syst,functions=( 'y0_t[5]', 'y0_s[5]'),
vars_dict={ 'r6':r6,'j6':j6,'dn1N2_dL': dn1N2_dL, 'dn2N2_dL':dn2N2_dL})
bvp_int.load_bcs({'a': {'1': y0_t[5]}, 'b': {'2': y0_s[5]}})
bvp_int.load_mesh(np.linspace(0, 0.15, 20))
sol = bvp_int.solve(tol=1e-9,verbose=2)

### Define a system Energy
syst = ['dT1_dL', '-dT2_dL'] # This is the dy/dt vector
bvp_int.define_system(syst,functions=( 'T0_t', 'T0_s'),
vars_dict={'U': U,'c': c,'sumReac': sumReac,'sumJicpi':sumJicpi,'sumFmcpitube':sumFmcpitube,
'sumFmcpisweep':sumFmcpisweep,'Tref':Tref, 'dT1_dL': dT1_dL, 'dT2_dL':dT2_dL})
bvp_int.load_bcs({'a': {'1': T0_t}, 'b': {'2': T0_s}})
bvp_int.load_mesh(np.linspace(0, 0.15, 20))
sol = bvp_int.solve(tol=1e-9,verbose=2)


Подробнее здесь: https://stackoverflow.com/questions/786 ... olving-bvp