两容器水重力流动模型溢流约束添加及仿真失败问题排查

问题背景 #

现有一个4变量4方程的两容器水重力流动Modelica模型，正常流动场景运行正常，但添加容器2满溢后流动停止、两容器压力均等的约束后，仿真在满溢时刻失败。

原模型代码 #

model TestOverflow
// two vessels are connected together with an opening on the bottom(e.g. pipe)
// the water shall flow by gravity from vessel 1 to vessel 2 until vessel 2 is fully filled
parameter Real HMax1(unit="m")=3; // height of vessel 1
parameter Real HMax2(unit="m")=1; // height of vessel 2
parameter Real A(unit="m2")=1; // cross section of vessels is equal to 1m2
constant Real rho(unit="kg/m3")=1000;
constant Real g(unit="m/s2")=9.81; // gravity
constant Real R=1.2; // flow resistance between two vessels

Real p1(unit="Pa"), p2(unit="Pa"); // pressure in vessel 1 and vessel 2
Real Q1(unit="kg/s"), Q2(unit="kg/s"); // mass flow into vessel 1,  vessel 2

initial equation
  p1=3*rho*g; // vessel 1 has 3 meter water level
  p2=0; // vessel 2 is empty
  
equation
Q1=-1/R * sqrt(2*rho*max(0,p1-p2)); // mass flow rate depends on the pressure difference and R=the flow resistance
// the flow shall stop when p1=p2 or less, so max(0,)

Q2=-Q1;  // conservation law
if p2<HMax2 * rho * g then  // normal flow
  Q1=A * der(p1) / g; // flow into vessel 1 increases the level = pressure/rho/g
  Q2=A * der(p2) / g;
else // vessel 2 is full, no flow, pressure p2 shall then become equal to p1!!! HOW TO REACH THIS???
  p2=p1;
  Q1=0;
end if;

end TestOverflow;

仿真错误信息 #

Simulation process failed. Exited with code 0xffffffffffffffff.
Error solving nonlinear system 8 at time 0.198323
nonlinear system 8 fails: at t=0.198323
[1] Real der(p2)(start=0, nominal=1)
Homotopy solver total pivot: Matrix (nearly) singular at time 0.198323.
under-determined linear system not solvable!
Homotopy solver total pivot: Matrix (nearly) singular at time 0.198323.
under-determined linear system not solvable!
Homotopy solver total pivot: Matrix (nearly) singular at time 0.198323.
under-determined linear system not solvable!
Error solving nonlinear system 8 at time 0.198323

问题根源 #

• 方程数量不匹配与系统欠定：在溢流分支（else块）中，原模型替换了两个微分方程，仅保留p2=p1和Q1=0，导致der(p1)和der(p2)失去约束，系统变量与有效方程数量不匹配，求解器无法处理欠定系统。
• 状态切换的导数不连续：直接替换微分方程会导致状态变量的导数在切换时刻突然无定义，引发求解器的奇异矩阵问题。
• 物理逻辑冲突：溢流时强制p2=p1和Q1=0与原流量方程重复约束，原流量方程在p1=p2时自然为0，额外约束反而引发矛盾。

修正方案 #

核心是保证两种状态下方程数量匹配，用物理合理的约束替代生硬分支，避免导数不连续。

修正后的模型代码 #

model TestOverflow_Fixed
// two vessels connected by bottom opening, flow stops when vessel 2 overflows and pressures equalize
parameter Real HMax1(unit="m")=3; // Height of vessel 1
parameter Real HMax2(unit="m")=1; // Height of vessel 2
parameter Real A(unit="m2")=1; // Cross-sectional area of vessels
constant Real rho(unit="kg/m3")=1000; // Water density
constant Real g(unit="m/s2")=9.81; // Gravitational acceleration
constant Real R=1.2; // Flow resistance between vessels

Real p1(unit="Pa"), p2(unit="Pa"); // Pressures in vessels 1 and 2
Real Q1(unit="kg/s"), Q2(unit="kg/s"); // Mass flows into vessels 1 and 2
Boolean isOverflow = p2 >= HMax2 * rho * g; // Auxiliary variable for overflow state

initial equation
  p1 = 3*rho*g; // Initial water level: 3m in vessel 1
  p2 = 0; // Vessel 2 starts empty

equation
  // Flow rate equation: naturally zero when p1 <= p2
  Q1 = -1/R * sqrt(2*rho*max(0, p1 - p2));
  Q2 = -Q1; // Mass conservation

  if not isOverflow then
    // Normal flow: pressure changes linked to mass flow
    Q1 = A * der(p1)/g;
    Q2 = A * der(p2)/g;
  else
    // Overflow state: vessel 2 pressure is capped at maximum
    p2 = HMax2 * rho * g;
    // Vessel 2 pressure stays constant, derivative is zero
    der(p2) = 0;
    // Vessel 1's pressure change is driven by outflow (which overflows vessel 2)
    Q1 = A * der(p1)/g;
  end if;

end TestOverflow_Fixed;

修正说明 #

• 方程数量匹配：两种状态下均保持4个有效方程，确保系统适定：
  • 正常流动：流量方程、质量守恒、两个微分方程
  • 溢流状态：流量方程、质量守恒、p2的约束方程、der(p2)=0的微分约束、Q1=A*der(p1)/g的微分方程
• 物理逻辑合理性：溢流时p2固定在最大允许压力，当p1高于该值时，水持续从容器1流出并溢出容器2，直到p1降至等于p2的最大值，流量自然停止，符合实际物理过程。
• 平滑状态切换：通过明确der(p2)=0替代移除微分方程，让求解器能连续处理状态切换，避免奇异矩阵问题。

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