我正在尝试模拟开始部分充满并逐渐填满的管道中的流动。一旦它满了，流量方程就需要改变。我已经尝试了以下代码，但是一旦管道已满，我就会得到除以零的错误。TestTrap1 是一个演示问题的演示。一旦管道已满，似乎某些方程式是不正确的。我正在使用 openmodelica 来运行它。
package FlowPackage
package Interfaces
extends Modelica.Icons.InterfacesPackage;
//Use this for a inlet outlet device
connector Fitting "Inlet port"
flow Modelica.SIunits.MassFlowRate q "Flow rate [Kg/s]";
Modelica.SIunits.Pressure P "Pressure";
annotation(defaultComponentName = "fitting", Icon(coordinateSystem(preserveAspectRatio = true, extent = {{-100, -100}, {100, 100}}), graphics = {Rectangle(extent = {{-100, 100}, {100, -100}}, lineColor = {0, 0, 255}, fillColor = {0, 0, 255}, fillPattern = FillPattern.Solid)}), Diagram(coordinateSystem(preserveAspectRatio = true, extent = {{-100, -100}, {100, 100}}), graphics = {Rectangle(extent = {{-40, 40}, {40, -40}}, lineColor = {0, 0, 255}, fillColor = {0, 0, 255}, fillPattern = FillPattern.Solid), Text(extent = {{-160, 110}, {40, 50}}, lineColor = {0, 0, 255}, textString = "%name")}));
end Fitting;
end Interfaces;
model SourceConstant
parameter Modelica.SIunits.MassFlowRate q = 0.015 "Flow rate [Kg/s]";
Interfaces.Fitting fitting annotation(Placement(visible = true, transformation(origin = {100, -6}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {100, -6}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
equation
-fitting.q = q;
annotation(Icon(graphics = {Rectangle(origin = {-3, 30}, extent = {{-93, 60}, {93, -60}}), Rectangle(origin = {0, 15}, fillColor = {0, 0, 255}, fillPattern = FillPattern.Solid, extent = {{-90, 45}, {90, -45}})}), Diagram);
end SourceConstant;
model ZeroP
parameter Modelica.SIunits.Pressure p = 0.00 "Pressure";
Interfaces.Fitting fitting annotation(Placement(visible = true, transformation(origin = {100, -6}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {100, -6}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
equation
fitting.P = p;
annotation(Icon(graphics = {Rectangle(origin = {-1, 2}, extent = {{-95, 94}, {95, -94}})}), Diagram);
end ZeroP;
model TestTrap1
TrapTemp trap1 ( level(start = 0.0509)) annotation(Placement(visible = true, transformation(origin = {10, 8}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
SourceConstant sourceConstant1(q = 0.03) annotation(Placement(visible = true, transformation(origin = {-50, 26}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
ZeroP zeroP1 (p = 101000) annotation(Placement(visible = true, transformation(origin = {64, 8}, extent = {{10, -10}, {-10, 10}}, rotation = 0)));
equation
connect(sourceConstant1.fitting, trap1.fitting_in) annotation(Line(points = {{-40, 26}, {-18, 26}, {-18, 8}, {0, 8}, {0, 8}}, color = {0, 0, 255}));
connect(trap1.fitting_out, zeroP1.fitting) annotation(Line(points = {{20, 8}, {54, 8}, {54, 8}, {54, 8}}, color = {0, 0, 255}));
annotation(Icon, Diagram, experiment(StartTime = 0, StopTime = 1450, Tolerance = 1e-6, Interval = 2.9));
end TestTrap1;
model TrapTemp
//Calculate the velocity for pipe segments
Interfaces.Fitting fitting_in annotation(Placement(visible = true, transformation(origin = {-100, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {-100, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
Interfaces.Fitting fitting_out annotation(Placement(visible = true, transformation(origin = {98, -2}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {98, -2}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
import m = Modelica.Math;
parameter Modelica.SIunits.Length L = 25000 "pipe length (m)";
parameter Modelica.SIunits.Radius R = 0.10194 / 2.0 "pipe inner radius (m)";
parameter Modelica.SIunits.Height roughness = 6.0e-5 "Average height of surface asperities (default: smooth steel pipe)";
/*Values below are liquid at 300K*/
parameter Modelica.SIunits.DynamicViscosity mu = 0.092977 "Dynamic Viscocity (cp)";
parameter Modelica.SIunits.Density row = 501.22 "Density (kg/m^3)";
Modelica.SIunits.Area A;
Modelica.SIunits.Length s(start = 0.158) "wetted perimeter";
Modelica.SIunits.Length rh "hydraulic radius";
Modelica.SIunits.Diameter De(start = 0.1) "eqivalent diameter [m]";
Modelica.SIunits.ReynoldsNumber Re(start = 3000);
Modelica.SIunits.Velocity v(start = 0.0001);
Modelica.SIunits.Length hf "Friction head loss";
Modelica.SIunits.Acceleration g = Modelica.Constants.g_n;
Modelica.SIunits.Acceleration a "Acceleration of flow";
Modelica.SIunits.CoefficientOfFriction f;
Real relative_roughness;
Real percent_full "Fraction of pipe is full [%]";
Modelica.SIunits.Mass mass;
Real percent_full "Fraction of pipe is full [%]";
Modelica.SIunits.Height level "Liquid level (m)";
equation
percent_full = level / R * 100.0;
relative_roughness = roughness / De;
a = der(v);
s = 2 * level + R;
rh = A / s;
De = 4 * rh;
//Protect against negative Re:
Re = De * abs(v) * row / mu;
//Head loss due to friction. Corrected for direction
hf = sign(v) * f * L * v ^ 2 / (2 * De * g);
f = 64 / Re ;
mass = A * row * L;
fitting_in.q + fitting_out.q = der(mass);
fitting_out.q = -v * A * row;
mass * a = A * (fitting_in.P - fitting_out.P - hf * row * g ) ;
A = R * level;
if noEvent(level >= R) then
//full pipe
fitting_in.q = -fitting_out.q;
else
//partially full pipe
fitting_in.P = fitting_out.P;
end if;
annotation(Icon, Diagram);
end TrapTemp;
annotation(Icon, Diagram);
end FlowPackage;
任何关于为什么这不起作用或如何解决它的建议将不胜感激。