使用 GEKKO 求解由于初始条件而具有共线方程的 ODE 系统

Question

我正在尝试求解以下代表圆锥形喷动床反应器的 ODE 系统，其具有 6 个未知数，即 4 个速度（us、ua、vs、va）、压力 p 和空隙率 eps_s。

我在 z=0 处的初始条件是 us = 24 m/s、ua = vs = va = 0 m/s、p = 360000 Pa 和 eps_s = 1。 然而，由于我的初始条件，系统的最后两个方程在 z=0 处共线，并且到目前为止我尝试过的所有求解器（包括 MATLAB 和 Julia 的显式和隐式求解器）都无法求解该系统。这就是为什么我想尝试看看是否可以使用 GEKKO 来解决我的系统问题。

不幸的是，GEKKO 也失败了，我想知道这是系统本身的原因还是我在实现过程中犯了错误。

import numpy as np
from gekko import GEKKO

m = GEKKO(remote=False)
m.time = np.linspace(0., 0.2,10)

# parameters
D0 = m.Param(0.06)
D_s = m.Param(0.04)
alpha_c = m.Param(1.152)
eps_a = m.Param(0.45)
eta_g = m.Param(1e-5)
rho_g = m.Param(1.3)
spheri = m.Param(1.)
d_p = m.Param(0.003)
rho_s = m.Param(2360.)
g = m.Param(9.81)

# variables and initial conditions
u_s = m.Var(value=24.)
u_a = m.Var(value=0.)
v_s = m.Var(value=0.)
v_a = m.Var(value=0.)
p = m.Var(value=360000.)
eps_s = m.Var(value=1.)

# additional
z = m.Var(value=m.time)
D_c = D0 + 2.*z*m.tan(alpha_c*0.5)

Re_p_a = m.Intermediate(rho_g * m.abs(v_a-u_a) * d_p * eps_a / eta_g)
C_D_a = m.if3(-1*Re_p_a, 24./Re_p_a * (1. + 0.15*Re_p_a**0.687), 0)
beta_Ergun_a = 150. * (1.-eps_a)**2 * eta_g / (eps_a*(d_p*spheri)**2) + 1.75*rho_g*(1.-eps_a)*m.abs(u_a-v_a) / (d_p*spheri)
beta_Wen_a = 0.75 * C_D_a * rho_g*(1.-eps_a)*m.abs(u_a-v_a) / (d_p*spheri) * eps_a**-2.65
phi_g_a = m.atan(150*1.75*(0.2-(1.-eps_a))) / np.pi + 0.5
beta_a = (1.-phi_g_a)*beta_Ergun_a + phi_g_a*beta_Wen_a

Re_p_s = m.Intermediate(rho_g * m.abs(v_s-u_s) * d_p * eps_s / eta_g)
C_D_s = 24./Re_p_s * (1. + 0.15*Re_p_s**0.687)
beta_Ergun_s = 150. * (1.-eps_s)**2 * eta_g / (eps_s*(d_p*spheri)**2) + 1.75*rho_g*(1.-eps_s)*m.abs(u_s-v_s) / (d_p*spheri)
beta_Wen_s = 0.75 * C_D_s * rho_g*(1.-eps_s)*m.abs(u_s-v_s) / (d_p*spheri) * eps_s**-2.65
phi_g_s = m.atan(150*1.75*(0.2-(1.-eps_s))) / np.pi + 0.5
beta_s = (1.-phi_g_s)*beta_Ergun_s + phi_g_s*beta_Wen_s

# equations
m.Equation(-1.*(D_c**2 -D_s**2)*u_a.dt() - D_s**2 /(eps_a) * (eps_s*u_s.dt() + u_s *eps_s.dt()) - 4. *u_a * D_c *m.tan(alpha_c*0.5) == 0)
m.Equation(-1.*(D_c**2 -D_s**2)*v_a.dt() - D_s**2 /(1.-eps_a)* ((1.-eps_s)*v_s.dt() - v_s *eps_s.dt()) - 4. *v_a * D_c *m.tan(alpha_c*0.5) == 0)

m.Equation(-1.* u_s.dt() *u_s*eps_s*rho_g - p*eps_s.dt() - eps_s *p.dt() - beta_s*(u_s-v_s) - g*eps_s*(rho_s-rho_g) == 0)
m.Equation(-1.* v_s.dt() *v_s*(1.-eps_s)*rho_s + p*eps_s.dt() - (1.-eps_s)*p.dt() + beta_s*(u_s-v_s) - g*(1.-eps_s)*(rho_s-rho_g) == 0)

m.Equation(-2.*rho_g*u_a*(D_c**2 -D_s**2)*u_a.dt() - 4.*rho_g *u_a**2 *D_c*m.tan(alpha_c*0.5) + rho_g*u_a*D_s**2 /eps_a * (eps_s*u_s.dt() + u_s*eps_s.dt()) - (D_c**2 -D_s**2)*p.dt() - 4.*p*D_c*m.tan(alpha_c*0.5) - (D_c**2 -D_s**2)*beta_a*(u_a-v_a) /eps_a - (D_c**2 -D_s**2)*g*(rho_s-rho_g) == 0)
m.Equation(-2.*rho_s*v_a*(D_c**2 -D_s**2)*v_a.dt() - 4.*rho_s *v_a**2 *D_c*m.tan(alpha_c*0.5) - rho_s*v_a*D_s**2 /(1.-eps_a) * ((1.-eps_s)*v_s.dt() - v_s*eps_s.dt()) - (D_c**2 -D_s**2)*p.dt() - 4.*p*D_c*m.tan(alpha_c*0.5) + (D_c**2 -D_s**2)*beta_a*(u_a-v_a) /(1.-eps_a) - (D_c**2 -D_s**2)*g*(rho_s-rho_g) == 0)

# solve
m.options.IMODE = 4
m.solve()
m.cleanup()

提前致谢。

Answer 1

如果删除

m.if3()

语句，则方程组求解成功。通过

m.if3()

语句，求解器确定该方程组不可行。

#C_D_a = m.if3(-1*Re_p_a, 24./Re_p_a * (1. + 0.15*Re_p_a**0.687), 0)
C_D_a = 24./Re_p_a * (1. + 0.15*Re_p_a**0.687)

我建议检查

m.if3()

条件，以确保

C_D_a

如果

Re_p_a<0

应使用零值。

Re_p_c

雷诺数计算没有相同的条件。最终解决方案中

Re_p_a

的值也显得非常小，但为正值，因此您可能不需要

m.if3()

条件。

Re_p_a = [0.0, 1.9048571431e-11, 1.4395918458e-11, 2.9893167235e-12, 
          8.6407817168e-13, 9.5340975304e-13, 9.020616017e-13, 
          4.5442810476e-12, 3.2721611684e-12, 1.4398004515e-12]

对于如此小的数字，

m.options.RTOL

（默认=

1e-6

）求解器收敛容差也可能存在问题。这是一个成功解决的完整脚本：

import numpy as np
from gekko import GEKKO

m = GEKKO(remote=False)
m.time = np.linspace(0., 0.2,10)

# parameters
D0 = m.Param(0.06)
D_s = m.Param(0.04)
alpha_c = m.Param(1.152)
eps_a = m.Param(0.45)
eta_g = m.Param(1e-5)
rho_g = m.Param(1.3)
spheri = m.Param(1.)
d_p = m.Param(0.003)
rho_s = m.Param(2360.)
g = m.Param(9.81)

# variables and initial conditions
u_s = m.Var(value=24.)
u_a = m.Var(value=0.)
v_s = m.Var(value=0.)
v_a = m.Var(value=0.)
p = m.Var(value=360000.)
eps_s = m.Var(value=1.)

# additional
z = m.Var(value=m.time)
D_c = D0 + 2.*z*m.tan(alpha_c*0.5)

Re_p_a = m.Intermediate(rho_g * m.abs(v_a-u_a) * d_p * eps_a / eta_g)
#C_D_a = m.if3(-Re_p_a+1e-10, 0, 24./Re_p_a * (1. + 0.15*Re_p_a**0.687))
C_D_a = 24./Re_p_a * (1. + 0.15*Re_p_a**0.687)
beta_Ergun_a = 150. * (1.-eps_a)**2 * eta_g / (eps_a*(d_p*spheri)**2) + 1.75*rho_g*(1.-eps_a)*m.abs(u_a-v_a) / (d_p*spheri)
beta_Wen_a = 0.75 * C_D_a * rho_g*(1.-eps_a)*m.abs(u_a-v_a) / (d_p*spheri) * eps_a**-2.65
phi_g_a = m.atan(150*1.75*(0.2-(1.-eps_a))) / np.pi + 0.5
beta_a = (1.-phi_g_a)*beta_Ergun_a + phi_g_a*beta_Wen_a

Re_p_s = m.Intermediate(rho_g * m.abs(v_s-u_s) * d_p * eps_s / eta_g)
C_D_s = 24./Re_p_s * (1. + 0.15*Re_p_s**0.687)
beta_Ergun_s = 150. * (1.-eps_s)**2 * eta_g / (eps_s*(d_p*spheri)**2) + 1.75*rho_g*(1.-eps_s)*m.abs(u_s-v_s) / (d_p*spheri)
beta_Wen_s = 0.75 * C_D_s * rho_g*(1.-eps_s)*m.abs(u_s-v_s) / (d_p*spheri) * eps_s**-2.65
phi_g_s = m.atan(150*1.75*(0.2-(1.-eps_s))) / np.pi + 0.5
beta_s = (1.-phi_g_s)*beta_Ergun_s + phi_g_s*beta_Wen_s

# equations
m.Equation(-1.*(D_c**2 -D_s**2)*u_a.dt() - D_s**2 /(eps_a) * (eps_s*u_s.dt() + u_s *eps_s.dt()) - 4. *u_a * D_c *m.tan(alpha_c*0.5) == 0)
m.Equation(-1.*(D_c**2 -D_s**2)*v_a.dt() - D_s**2 /(1.-eps_a)* ((1.-eps_s)*v_s.dt() - v_s *eps_s.dt()) - 4. *v_a * D_c *m.tan(alpha_c*0.5) == 0)

m.Equation(-1.* u_s.dt() *u_s*eps_s*rho_g - p*eps_s.dt() - eps_s *p.dt() - beta_s*(u_s-v_s) - g*eps_s*(rho_s-rho_g) == 0)
m.Equation(-1.* v_s.dt() *v_s*(1.-eps_s)*rho_s + p*eps_s.dt() - (1.-eps_s)*p.dt() + beta_s*(u_s-v_s) - g*(1.-eps_s)*(rho_s-rho_g) == 0)

m.Equation(-2.*rho_g*u_a*(D_c**2 -D_s**2)*u_a.dt() - 4.*rho_g *u_a**2 *D_c*m.tan(alpha_c*0.5) + rho_g*u_a*D_s**2 /eps_a * (eps_s*u_s.dt() + u_s*eps_s.dt()) - (D_c**2 -D_s**2)*p.dt() - 4.*p*D_c*m.tan(alpha_c*0.5) - (D_c**2 -D_s**2)*beta_a*(u_a-v_a) /eps_a - (D_c**2 -D_s**2)*g*(rho_s-rho_g) == 0)
m.Equation(-2.*rho_s*v_a*(D_c**2 -D_s**2)*v_a.dt() - 4.*rho_s *v_a**2 *D_c*m.tan(alpha_c*0.5) - rho_s*v_a*D_s**2 /(1.-eps_a) * ((1.-eps_s)*v_s.dt() - v_s*eps_s.dt()) - (D_c**2 -D_s**2)*p.dt() - 4.*p*D_c*m.tan(alpha_c*0.5) + (D_c**2 -D_s**2)*beta_a*(u_a-v_a) /(1.-eps_a) - (D_c**2 -D_s**2)*g*(rho_s-rho_g) == 0)

# solve
m.options.IMODE = 6
m.solve()
m.cleanup()

另一件事需要注意的是

m.abs()

函数对于基于梯度的求解器来说不是连续可微的。使用

m.abs3()

可能有助于找到解决方案。与任何应用程序一样，最好使用示例时间点来验证您的预期解决方案。