Modelica: Program a static solver to converge to a solution during the simulation. CFD example

Question

Generic question:

I work on Modelica with OpenModelica. I would like to program my own solver in Modelica to converge to a solution because I think OpenModelica can't solve my problem. Is it possible?

---

Specific case:

I developed a model of a separation of air in fluid dynamic component. My model is based on a data table which provide the coefficient of pressure loss for each branche depending of the velocities in the component and the parameter "section". My model works well when I use MassFlow sources but not when I use only Pressure source. See the pictures below to understand the connection with my component:

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

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

The code of my model "separation" is this one:

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;

It's a "nonlinear system" error when I connect this model to 3 Pressure components (set to Patm +10000Pa for the source and Patm for the sink). This model works well with the MassFlow sink. Why? Should I develop my own solver to solve it? If yes how?

Answer 1

Your question is not completely clear. But, I think what you are asking if given a value of dp, at any time, can modelica solve for v, when there is a table for the relationship of K and v? If I understand your variables correctly, this should do it:

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;

When I run this, it solves for v, given a value of dp. Is that what you want?

Answer 2

OpenModelica gives a hint that the original model is bad:

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 is defined in both an initial algorithm and an algorithm section. initial algorithms are only supposed to be used for parameters or states (both der(x) and x need to be specified in the initial step). Luckily, you can specify v(start=0), and v will be initialized to 0 at the start of an algorithm section that assigns to v.

After that is resolved, the model compiles and simulates. Although as others stated, algorithm sections are ugly and should be avoided at all costs in Modelica.

Answer 3

I'm not sure I understand the problem 100%. You said that you could not express it analytically. But can you express the relationship between K and v as a Modelica function? If so, you can use the inverse annotation to supply an inverse function. In this way, the tool will use the most efficient version of the function and can avoid having to do any non-linear iteration. But this assumes you can formulate the inverse function. It wasn't clear to me, from your explanation, whether you could or not.