基于MSL TemplateMedium的自定义常密度介质初始化错误排查

问题简述 #

基于MSL媒体包TemplateMedium创建的自定义介质，原本可通过MSL介质测试，但在带阀门的模型中报错。改为常密度介质后，连基础测试模型都无法通过初始化。推测原介质通过密度与压力的关联抵消管道压力增量，维持压力等于p_start；而常密度介质缺失该机制，引发初始化矛盾。

报错信息 #

The initialization problem is inconsistent due to the following equation:
0 != -10000 = volume.p_start - volume.medium.p

相关Modelica代码 #

package MyPackageSO

  package SimpleMedia "Simple media with constant density and linear enthalpy"
    extends Modelica.Media.Interfaces.PartialMedium(final mediumName = "Diesel", final substanceNames = {mediumName}, final singleState = false, final reducedX = true, final fixedX = true, Temperature(min = 273, max = 373, start = 298), p_default = 1e5, reference_p = 1e5, reference_T = 298, reference_X = {1}, AbsolutePressure(min = 1e-6, max = 3e8, start = 1e5),ThermoStates = Modelica.Media.Interfaces.Choices.IndependentVariables.pT);
      // Provide medium constants here
      constant SpecificHeatCapacity cp_const = 1350 "Constant specific heat capacity at constant pressure -- using approx value at 300K";
  
    redeclare model extends BaseProperties(final standardOrderComponents = true) "Base properties of medium"
      equation
        d = density(state);
        h = specificEnthalpy(state);
        u = h - p/d;
        MM = 0.025;
        R_s = 0;
        state.p = p;
        state.T = T;
    end BaseProperties;
  
    redeclare replaceable record ThermodynamicState "A selection of variables that uniquely defines the thermodynamic state"
      extends Modelica.Icons.Record;
      AbsolutePressure p "Absolute pressure of medium";
      Temperature T "Temperature of medium";
      annotation(
        Documentation(info = "<html>
    
    </html>"));
    end ThermodynamicState;
  
    redeclare function extends setState_pTX "Return thermodynamic state as function of p, T and composition X or Xi"
      extends Modelica.Icons.Function;
      algorithm
        state := ThermodynamicState(p = p, T = T);
      annotation(
        Inline = true);
    end setState_pTX;
    
    redeclare function extends setState_phX "Return thermodynamic state as function of p, h and composition X or Xi"
      extends Modelica.Icons.Function;
      algorithm
        state := ThermodynamicState(p = p, T = h / cp_const);
      annotation(
        Inline = true);
    end setState_phX;
      
    redeclare function extends pressure "Return the Pressure"
    algorithm
      p := state.p;
    end pressure;
    
    redeclare function extends temperature "Return temperature"
      input ThermodynamicState state "Thermodynamic state record";
      output Temperature T "Temperature";
    algorithm
      T := state.T;
    end temperature;
    
    redeclare function extends density "Return the density"
    algorithm
      d := 871;
    end density;
    
    redeclare function extends specificEnthalpy "Return specific enthalpy"
      extends Modelica.Icons.Function;
      input ThermodynamicState state "Thermodynamic state record";
      output SpecificEnthalpy h "Specific enthalpy";
      algorithm  
        h := cp_const * (state.T - reference_T);
      annotation(Documentation(info = "<html></html>"));
    end specificEnthalpy;      
    
        redeclare function extends specificInternalEnergy
      "Return specific internal energy"
      extends Modelica.Icons.Function;
    algorithm
      u := cp_const*(state.T - reference_T);
      annotation (Documentation(info="<html>
  <p>
  This function computes the specific internal energy of the fluid, but neglects the (small) influence of the pressure term p/d.
  </p>
  </html>"));
    end specificInternalEnergy;
    
    redeclare function extends dynamicViscosity "Return dynamic viscosity"
      algorithm
        eta := 0.014773467;
      annotation(
        Documentation(info = "<html>
    
    </html>"));
    end dynamicViscosity;
    
    redeclare function extends thermalConductivity "Return thermal conductivity"
      algorithm
        lambda := 0.1;
      annotation(
        Documentation(info = "<html>
    
    </html>"));
    end thermalConductivity;
    
    redeclare function extends specificHeatCapacityCp "Return specific heat capacity at constant pressure"
      algorithm
        cp := 1350;
      annotation(
        Documentation(info = "<html>
    
    </html>"));
    end specificHeatCapacityCp;
    
    redeclare function extends isentropicExponent "Return isentropic exponent"
        extends Modelica.Icons.Function;
    
      algorithm
        gamma := 1.16;
      annotation(
        Documentation(info = "<html>
    
    </html>"));
    end isentropicExponent;
    
    redeclare function extends velocityOfSound "Return velocity of sound"
        extends Modelica.Icons.Function;
      algorithm
        a := 1100;
      annotation(
        Documentation(info = "<html>
    
    </html>"));
    end velocityOfSound;
  
  end SimpleMedia;

  model TestOfMyMedium
    extends Modelica.Media.Examples.Utilities.PartialTestModel(redeclare package Medium = MyPackage.SimpleMedia, volume(p_start = 1e8), ambient(p_ambient = 1e8));
  end TestOfMyMedium;
end MyPackageSO;

修改方案 #

1. 修正setState_phX函数的温度计算 #

原函数通过焓计算温度时忽略了参考温度reference_T，与specificEnthalpy的定义逻辑断裂。根据h = cp_const*(T - reference_T)反推，温度应为T = reference_T + h/cp_const，修改后代码：

redeclare function extends setState_phX "Return thermodynamic state as function of p, h and composition X or Xi"
  extends Modelica.Icons.Function;
  algorithm
    state := ThermodynamicState(p = p, T = reference_T + h / cp_const);
  annotation(
    Inline = true);
end setState_phX;

2. 对齐specificInternalEnergy与BaseProperties的方程 #

BaseProperties中已定义u = h - p/d，但原specificInternalEnergy函数忽略了p/d项，导致方程冲突。修改函数使其与BaseProperties的定义一致：

redeclare function extends specificInternalEnergy
  "Return specific internal energy"
  extends Modelica.Icons.Function;
algorithm
  u := specificEnthalpy(state) - pressure(state)/density(state);
  annotation (Documentation(info="<html>
<p>
Specific internal energy calculated as h - p/d, consistent with BaseProperties equation.
</p>
</html>"));
end specificInternalEnergy;

或直接展开计算（效率更高）：

redeclare function extends specificInternalEnergy
  "Return specific internal energy"
  extends Modelica.Icons.Function;
algorithm
  u := cp_const*(state.T - reference_T) - state.p/871;
  annotation (Documentation(info="<html>
<p>
Specific internal energy calculated as h - p/d, consistent with BaseProperties equation.
</p>
</html>"));
end specificInternalEnergy;

3. 统一初始化参数（可选优化） #

测试模型中p_start设为1e8，而介质的p_default为1e5，建议将介质的p_default与测试模型的p_start保持一致，减少初始化歧义：

extends Modelica.Media.Interfaces.PartialMedium(..., p_default = 1e8, ...);

原理说明 #

初始化矛盾的核心是介质状态函数实现不一致：setState_phX的温度计算错误导致焓与温度映射断裂，specificInternalEnergy与BaseProperties方程冲突，使得求解器无法找到满足所有方程的初始值。修正后，所有热力学函数逻辑统一，求解器可正常完成初始化。

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