使用GEKKO求解因初始条件导致方程共线的ODE系统

我正在求解一组描述锥形喷动床反应器的ODE系统，包含6个未知量：

• 4个速度：$u_s$（中心区气相速度）、$u_a$（环形区气相速度）、$v_s$（中心区固相速度）、$v_a$（环形区固相速度）
• 压力$p$
• 空隙率$\varepsilon_s$

初始条件为$z=0$时：

• $u_s=24\ \text{m/s}$，$u_a=v_s=v_a=0\ \text{m/s}$
• $p=360000\ \text{Pa}$，$\varepsilon_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()

问题排查方向 #

1. 初始条件奇异性
   当$\varepsilon_s=1$时，中心区无固相，导致$(1-\varepsilon_s)=0$，固气相互作用项消失，引发数值奇异。物理上喷动床初始状态不会是完全气相，建议将$\varepsilon_s$初始值设为略小于1（如0.999），避开奇异点。

2. 轴向坐标定义错误
   代码中z = m.Var(value=m.time)的逻辑错误：$z$是轴向坐标，需通过运动学方程与速度关联（如$\frac{dz}{dt}=u_s$），而非直接绑定时间轴，否则会破坏系统物理一致性。

3. 曳力模型切换逻辑错误
   C_D_a的if3触发条件为-1*Re_p_a，但雷诺数用绝对值计算，不可能为负，导致C_D_a始终为0，环形区曳力模型仅保留Ergun项，与实际物理模型不符。应改为按临界雷诺数（如$Re<1$）切换模型。

4. 方程线性相关性
   初始条件下$u_a=v_a=0$，最后两个环形区动量方程大部分项为0，仅剩压力梯度和重力项，可能导致方程线性相关，出现共线。需重新验证方程推导的正确性，排查冗余或错误项。

内容的提问来源于stack exchange，提问作者MoDi