OpenModelica Fluid库中流动变量m_flow的导数求解问题

问题原因 #

Modelica Fluid库中，端口的m_flow（质量流量）默认是代数变量，而非状态变量。代数变量由系统约束方程实时计算得出，没有直接的微分方程定义其变化率。直接对代数变量求导der(m_flow)会破坏微分代数方程(DAE)的求解结构，导致求解器无法处理无定义的导数项，进而触发报错。

可行解决方案 #

方案1：使用Derivative块近似导数（推荐） #

利用Modelica标准库中的Modelica.Blocks.Continuous.Derivative块，将代数变量的导数转换为状态变量跟踪，让求解器可以合法处理。修改测试模型如下：

model Fluid_test
  replaceable package Medium = Modelica.Media.Water.StandardWaterOnePhase
    constrainedby Modelica.Media.Interfaces.PartialMedium;
    
  Modelica.Fluid.Sources.FixedBoundary boundary(
      nPorts = 1,
      use_T=true,
      T=Modelica.Units.Conversions.from_degC(20),
      p=600000,
      redeclare package Medium = Medium) annotation(
      Placement(transformation(origin = {-56, -18}, extent = {{-10, -10}, {10, 10}})));
      
  Modelica.Fluid.Sources.FixedBoundary boundary1(
      p=600000,
      T=400,
      nPorts=2,
      redeclare package Medium = Medium)
      annotation (Placement(transformation(origin = {-18, 44}, extent = {{100, -40}, {80, -20}})));

  Modelica.Fluid.Machines.PrescribedPump pump(
      checkValve=true,
      checkValveHomotopy = Modelica.Fluid.Types.CheckValveHomotopyType.Closed,
      N_nominal=1200,
      redeclare function flowCharacteristic =
          Modelica.Fluid.Machines.BaseClasses.PumpCharacteristics.quadraticFlow (
            V_flow_nominal={0,0.25,0.5}, head_nominal={100,60,0}),
      use_N_in=true,
      nParallel=1,
      energyDynamics=Modelica.Fluid.Types.Dynamics.FixedInitial,
      V(displayUnit="l") = 0.05,
      massDynamics=Modelica.Fluid.Types.Dynamics.FixedInitial,
      redeclare package Medium = Medium,
      p_b_start=600000,
      T_start=400) annotation(
      Placement(transformation(origin = {-10, 14}, extent = {{-10, -10}, {10, 10}})));
      
  Modelica.Blocks.Sources.Constant const annotation(
      Placement(transformation(origin = {-8, 58}, extent = {{-10, -10}, {10, 10}})));
  
  // 添加Derivative块来跟踪m_flow的导数
  Modelica.Blocks.Continuous.Derivative der_m_flow(
      k=1, T=0.01) annotation(
      Placement(transformation(origin = {20, 14}, extent = {{-10, -10}, {10, 10}})));
  Real TestValue;

equation
  TestValue = der_m_flow.y;
  connect(boundary.ports[1], pump.port_a) annotation(
      Line(points = {{-46, -18}, {-20, -18}, {-20, 14}}, color = {0, 127, 255}));
  connect(pump.port_b, boundary1.ports[1]) annotation(
      Line(points = {{0, 14}, {62, 14}}, color = {0, 127, 255}));
  connect(const.y, pump.N_in) annotation(
      Line(points = {{4, 58}, {-10, 58}, {-10, 24}}, color = {0, 0, 127}));
  // 将m_flow连接到Derivative块的输入
  connect(pump.port_a.m_flow, der_m_flow.u) annotation(
      Line(points = {{0, 14}, {10, 14}}, color = {0, 0, 127}));
  annotation(
      Diagram);
end Fluid_test;

注：Derivative块的T参数可调整滤波时间常数，平衡导数的平滑性和响应速度。

方案2：添加缓冲容积间接推导导数 #

在泵的端口前串联一个极小的固定容积（Modelica.Fluid.Volumes.FixedVolume），利用容积的质量变化率与流入流量的关系，通过容积质量的二阶导数得到流量的导数：

Modelica.Fluid.Volumes.FixedVolume small_volume(
    V=1e-6, // 极小容积，避免影响系统特性
    redeclare package Medium=Medium,
    p_start=600000,
    T_start=Modelica.Units.Conversions.from_degC(20)) annotation(...);

// 连接关系：boundary -> small_volume -> pump.port_a
equation
// small_volume的质量变化率等于流入的m_flow
// der(small_volume.m) = pump.port_a.m_flow
// 因此流量的导数等于质量的二阶导数
TestValue = der(der(small_volume.m));

这种方法更贴合Fluid库的物理逻辑，但需要注意容积大小的设置，避免引入不必要的动态延迟。

方案3：数值差分近似（仿真阶段用） #

如果仅在仿真后处理时需要导数，可采用数值差分计算，比如在结果文件中用相邻时刻的流量差除以时间步长：

TestValue = (m_flow[t] - m_flow[t-Δt])/Δt

这种方法无需修改模型，但精度受时间步长影响，适合离线分析场景。

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