Modelica:对静态求解器进行编程，使其在模拟过程中收敛到一个解.CFD的例子

Question

通用性问题:

我使用OpenModelica开发Modelica。我想在Modelica编程我自己的求解器收敛到一个解决方案，因为我认为OpenModelica不能解决我的问题。这可能吗?

---

具体的案例:

我建立了一个流体动力成分中空气分离的模型。我的模型是基于一个数据表，该数据表提供了每个分支的压力损失系数，这取决于组件中的速度和参数"截面"。我的模型工作得很好，当我使用MassFlow源，但不是当我只使用压力源。查看下面的图片来理解与我的组件的连接:

http://www.casimages.com/img.php?i=140620024048704083.png

http://www.casimages.com/img.php?i=140620024137384886.png

我的模型"separation"的代码是这样的:

model separation
 replaceable package Medium = Modelica.Media.Interfaces.PartialMedium "Medium in the component";
 Modelica.Fluid.Interfaces.FluidPort_a port_a(redeclare package Medium = Medium);
 Modelica.Fluid.Interfaces.FluidPort_b port_b2(redeclare package Medium = Medium);
 Modelica.Fluid.Interfaces.FluidPort_b port_b1(redeclare package Medium = Medium);
 Modelica.Blocks.Tables.CombiTable2D coeff_PDC1(table = [0, 0, 0.4, 0.5, 0.6, 0.7, 0.8, 1; 0, 1, 1, 1, 1, 1, 1, 1; 0.1, 0.81, 0.81, 0.81, 0.81, 0.81, 0.81, 0.81; 0.2, 0.64, 0.64, 0.64, 0.64, 0.64, 0.64, 0.64; 0.3, 0.5, 0.5, 0.52, 0.52, 0.5, 0.5, 0.5; 0.4, 0.36, 0.36, 0.4, 0.38, 0.37, 0.36, 0.36; 0.5, 0.25, 0.25, 0.3, 0.28, 0.26, 0.25, 0.25; 0.6, 0.16, 0.16, 0.23, 0.2, 0.18, 0.16, 0.16; 0.8, 0.04, 0.04, 0.16, 0.12, 0.07, 0.04, 0.04; 1, 0.001, 0.001, 0.2, 0.1, 0.05, 0.001, 0.001; 1.2, 0.07, 0.07, 0.36, 0.21, 0.14, 0.07, 0.07; 1.4, 0.39, 0.39, 0.78, 0.59, 0.49, 50, 50; 1.6, 0.9, 0.9, 1.36, 1.15, 50, 50, 50; 1.8, 1.78, 1.78, 2.43, 50, 50, 50, 50; 2, 3.2, 3.2, 4, 50, 50, 50, 50]);
 Modelica.Blocks.Tables.CombiTable1Ds coeff_PDC2( table = [0.1, 1; 0.2, 1; 0.3, 1; 0.4, 1; 0.5, 1; 0.6, 1; 0.8, 1; 1,1; 1.2, 1; 1.4, 1; 1.6, 1; 1.8, 1; 2, 1]);
 Real dp_b1 "en Pa, perte de charge entre les ports a1 et b";
 Real dp_b2 "en Pa, perte de charge entre les ports a2 et b";
 Modelica.SIunits.Velocity v_a(start = 0);
 Modelica.SIunits.Velocity v_b1(start = 0);
 Modelica.SIunits.Velocity v_b2(start = 0);
 parameter Real rho=1.2;
 parameter Modelica.SIunits.Area surface_b1 = 1;
 parameter Modelica.SIunits.Area surface_b2 = 1;
 parameter Modelica.SIunits.Area surface_a = 2;
equation 
 coeff_PDC1.u1 = if noEvent(abs(v_a) > 0) then v_b1/v_a else 1;
 coeff_PDC1.u2 = surface_b1/surface_a;
 coeff_PDC2.u = if noEvent(abs(v_a) > 0) then v_b2/v_a else 1;
 v_a = abs(port_a.m_flow)/rho/surface_a;
 v_b1 = abs(port_b1.m_flow)/rho/surface_b1;
 v_b2 = abs(port_b2.m_flow)/rho/surface_b2;
 port_b1.p - port_a.p = dp_b1;
 dp_b1 = 1/2*coeff_PDC1.y*port_b1.m_flow^2/surface_b1^2/rho;
 port_b2.p - port_a.p = dp_b2;
 dp_b2 = 1/2*coeff_PDC2.y[1]*port_b2.m_flow^2/surface_b2^2/rho;
 port_b1.m_flow + port_b2.m_flow + port_a.m_flow = 0;
 port_b1.Xi_outflow = inStream(port_a.Xi_outflow);
 port_b2.Xi_outflow = inStream(port_a.Xi_outflow);
 port_a.Xi_outflow = inStream(port_b1.Xi_outflow);
 port_b1.C_outflow = inStream(port_a.C_outflow);
 port_b2.C_outflow = inStream(port_a.C_outflow);
 port_a.C_outflow = inStream(port_b1.C_outflow);
 port_b1.h_outflow = inStream(port_a.h_outflow);
 port_b2.h_outflow = inStream(port_a.h_outflow);
 port_a.h_outflow = inStream(port_b1.h_outflow);
end separation;

当我将此模型连接到3个压力组件时，这是一个"非线性系统"错误(源设置为Patm +10000Pa，汇设置为Patm)。这个模型可以很好地与MassFlow接收器一起工作。为什么?我应该开发自己的求解器来解决它吗?如果是，怎么做?

Answer 1

你的问题不是很清楚。但是，我想你问的是如果给定一个dp值，在任何时候，当有K和v的关系表时，是否可以建模求解v?如果我对变量的理解正确的话，应该是这样的:

model test5
  Real v "fluid velocity";
  Real K "Pressure loss coefficient";
  Real dp =0.5 "Pressure drop"; //Example value
  Real rho = 1.0 "density" ;
  Modelica.Blocks.Tables.CombiTable2D coeff_PDC1(table = [0, 0, 0.4, 0.5, 0.6, 0.7, 0.8, 1; 0, 1, 1, 1, 1, 1, 1, 1; 0.1, 0.81, 0.81, 0.81, 0.81, 0.81, 0.81, 0.81; 0.2, 0.64, 0.64, 0.64, 0.64, 0.64, 0.64, 0.64; 0.3, 0.5, 0.5, 0.52, 0.52, 0.5, 0.5, 0.5; 0.4, 0.36, 0.36, 0.4, 0.38, 0.37, 0.36, 0.36; 0.5, 0.25, 0.25, 0.3, 0.28, 0.26, 0.25, 0.25; 0.6, 0.16, 0.16, 0.23, 0.2, 0.18, 0.16, 0.16; 0.8, 0.04, 0.04, 0.16, 0.12, 0.07, 0.04, 0.04; 1, 0, 0, 0.2, 0.1, 0.05, 0, 0; 1.2, 0.07, 0.07, 0.36, 0.21, 0.14, 0.07, 0.07; 1.4, 0.39, 0.39, 0.78, 0.59, 0.49, 50, 50; 1.6, 0.9, 0.9, 1.36, 1.15, 50, 50, 50; 1.8, 1.78, 1.78, 2.43, 50, 50, 50, 50; 2, 3.2, 3.2, 4, 50, 50, 50, 50]);
equation
  coeff_PDC1.u2 = 0.5;
  dp = 0.5 * K * rho * v ^ 2;
  v = coeff_PDC1.u1;
  K = coeff_PDC1.y;
end test5;

当我运行这个时，它解出v，给定值dp。这就是你想要的吗?

Answer 2

OpenModelica给出了原始模型不好的提示:

Error: Initialization problem is structurally singular, error found sorting equations 
  1: algorithm
    v := 0.0;
  2: algorithm
    K := 1.0;
    while v > 0.01 + $PRE.v loop
      v := (10.0 / K) ^ 0.5;
      K := 1.0 + v ^ 2.0;
    end while;

v在初始算法和算法章节中都有定义。初始算法仅用于参数或状态(der(x)和x都需要在初始步骤中指定)。幸运的是，您可以指定v(start=0)，并且v将在分配给v的算法部分开始时初始化为0。

解析后，模型进行编译和模拟。尽管正如其他人所说，算法节很难看，应该在Modelica中不惜一切代价避免使用。

Answer 3

我不确定我完全理解了这个问题。你说过你无法用分析的方式表达。但K和v之间的关系能否用Modelica函数来表示呢?如果是这样，您可以使用逆注释来提供逆函数。通过这种方式，该工具将使用函数的最有效版本，并且可以避免进行任何非线性迭代。但前提是你可以写出逆函数。从你的解释中，我不清楚你是否可以。