OpenModelica中DASSL求解PCM与外界换热的液相分数异常问题

一、DASSL求解器与时间导数的工作机制 #

DASSL是OpenModelica中针对**微分代数方程组（DAE）**的刚性隐式求解器，核心基于向后差分公式（BDF）实现：

• 它会先对输入的DAE进行预处理，自动完成降索引操作，把高阶微分方程转化为一阶残差形式；
• 求解过程中，通过迭代法同时处理微分方程的时间导数项和代数约束，无需手动拆分方程；
• 对于相变这类刚性系统（温度不变但潜热项导数剧烈变化），DASSL能通过隐式计算稳定捕捉系统状态，但对不连续条件的处理依赖事件触发机制，默认不会主动识别阶跃式的状态切换。

二、液相分数fl_xy异常的核心原因 #

1. 数值求解的精度误差：当温度接近熔点T_fusion时，求解器的截断误差可能导致计算温度略高于/低于熔点，但模型中阶跃式的液相分数逻辑（严格等于才允许0-1区间）无法适配这种微小波动；
2. 不连续条件未触发事件：如果模型中液相分数的切换仅靠条件判断而非事件触发，DASSL会将其视为连续函数处理，导致数值解偏离物理预期；
3. DAE约束一致性问题：空间离散后，温度与液相分数的耦合约束可能存在矛盾——比如当温度不在熔点时，未通过代数约束强制液相分数锁定为0或1，求解器可能给出中间值。

三、问题解决方案 #

1. 给相变逻辑添加数值鲁棒过渡 #

把严格阶跃的液相分数改为带微小过渡区间的连续函数，避免求解器处理不连续点时出错：

function LiquidFraction
  input Real T;
  input Real T_fusion;
  input Real deltaT = 0.1; // 可调整的过渡区间，平衡物理精度与数值稳定性
  output Real fl;
algorithm
  if T < T_fusion - deltaT/2 then
    fl := 0;
  elseif T > T_fusion + deltaT/2 then
    fl := 1;
  else
    fl := (T - (T_fusion - deltaT/2))/deltaT; // 线性过渡
  end if;
end LiquidFraction;

在Octadecanol模型中调用这个函数计算液相分数，既符合物理规律，又给求解器提供连续的计算区间。

2. 用事件触发强制约束一致性 #

在PCM_ex模型的方程部分，添加事件触发逻辑，确保温度偏离熔点时液相分数被强制锁定：

equation
  // 温度低于熔点时，强制液相分数为0
  when T_xy < T_fusion then
    fl_xy = 0;
  end when;
  // 温度高于熔点时，强制液相分数为1
  when T_xy > T_fusion then
    fl_xy = 1;
  end when;
  // 原有热平衡方程（保留你的换热与离散化逻辑）
  rho*cp*(1 - fl_xy)*der(T_xy) + rho*L*der(fl_xy) = ...;

when语句会触发求解器在温度穿越熔点时重新初始化状态，严格满足物理约束。

3. 调整DASSL求解器参数 #

在PCM_ex_test模型的实验配置中，开启事件检测并提高求解精度：

experiment(
  StopTime = 3600,
  Tolerance = 1e-6, // 减少数值误差对温度判断的影响
  Interval = 10,
  solver = "DASSL",
  solverOptions = "enableEventDetection=true; maxNonlinearIterations=50"
);

强制求解器识别温度穿越熔点的事件，避免跳过关键状态切换点。

4. 验证离散化方程的耦合逻辑 #

检查x、y方向离散后每个网格单元的热平衡方程：

• 确认潜热项rho*L*der(fl_xy)仅在温度等于熔点时生效；
• 非相变阶段，液相分数的导数der(fl_xy)需被约束为0，避免求解器给出无意义的中间值。

四、分步验证流程 #

1. 单独测试Octadecanol模型：给定点温度输入（低于、等于、高于熔点），检查液相分数输出是否符合预期；
2. 简化PCM_ex模型：去掉外界换热，只模拟PCM的纯相变过程，验证fl_xy随温度的变化逻辑；
3. 逐步添加外界换热与空间离散化，定位问题是否出在离散后的耦合方程中。

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