我试图使用Modelica对由弹性管道组成的系统进行建模。
这个模型应该更容易理解,因为它是一个平面模型,不从其他类扩展。这一点很重要,因为即使没有Modelica专有技术,我的同事也可以理解该模型,我可以说服他们,Modelica是实现我们目标的适当工具!
作为一个测试用例,我使用带有阶跃信号(水锤)的质量流源。
我的模型给出的结果与Modelica不同。流体成分。
这个
测试系统
看起来像这样:
11个单元格的结果
这是:
如您所见,MSL分量的压力峰值较高,频率/周期不同。当我选择更多的单元格时,误差就会变小。
我很确定我使用的是完全相同的方程。
这可能是数字原因(我尝试使用标称值)的原因吗?
我的管道型号代码
model Pipe_FVM_staggered
// Import
import SI = Modelica.SIunits;
import Modelica.Constants.pi;
// Medium
replaceable package Medium = Modelica.Media.Interfaces.PartialMedium "Medium in the component"
annotation (choicesAllMatching = true);
// Interfaces, Ports
Modelica.Fluid.Interfaces.FluidPort_a port_a(redeclare package Medium = Medium) annotation (Placement(transformation(extent={{-110,-10},{-90,10}})));
Modelica.Fluid.Interfaces.FluidPort_b port_b(redeclare package Medium = Medium) annotation (Placement(transformation(extent={{90,-10},{110,10}})));
// Parameters
parameter Integer n(min=2) = 3 "Number of cells"; // No effect yet, only for icon
parameter SI.Length L = 1 "Length";
parameter SI.Diameter D = 0.010 "Diameter";
parameter SI.Height R = 2.5e-5 "Roughness";
parameter Boolean use_fixed_zeta = false "Use fixed zeta value instead of Moody chart";
parameter SI.CoefficientOfFriction zeta = 1;
// Initialization
parameter Medium.Temperature T_start = 293.15 "Start temperature" annotation(Dialog(tab="Initialization"));
parameter Medium.MassFlowRate mflow_start = 1 "Start mass flow rate in design direction" annotation(Dialog(tab="Initialization"));
parameter Medium.AbsolutePressure p_a_start = 2e5 "Start pressure p[1] at design inflow" annotation(Dialog(tab="Initialization"));
parameter Medium.AbsolutePressure p_b_start = 1e5 "Start pressure for p[n+1] at design outflow" annotation(Dialog(tab="Initialization"));
// parameter Medium.AbsolutePressure p_start = (p_a_start + p_b_start)/2 annotation(Dialog(tab="Initialization"));
parameter Medium.AbsolutePressure p_start[:] = linspace(p_a_start, p_b_start, n) annotation(Dialog(tab="Initialization"));
// parameter Medium.SpecificEnthalpy h_start[:] = Medium.specificEnthalpy_pTX(p_start, T_start, Medium.X_default);
parameter Medium.SpecificEnthalpy h_start = Medium.specificEnthalpy_pTX((p_a_start + p_b_start)/2, T_start, Medium.X_default) annotation(Dialog(tab="Initialization"));
parameter SI.AbsolutePressure dp_nominal = 1e5;
parameter SI.MassFlowRate m_flow_nominal = 1;
// Variables general
SI.Length dL = L/n;
SI.Area A(nominal=0.001) = D^2*pi/4;
SI.Volume V = A * dL;
// Variables cell centers: positiv in direction a -> b
Medium.AbsolutePressure p[n](start = p_start, each stateSelect=StateSelect.prefer) annotation(Dialog(tab="Initialization", showStartAttribute=true, enable=false));
Medium.SpecificEnthalpy h[n](each start = h_start, each stateSelect=StateSelect.prefer) annotation(Dialog(tab="Initialization", showStartAttribute=true, enable=false));
Medium.ThermodynamicState state[n] = Medium.setState_phX(p,h);
SI.Mass m[n] = rho .* V;
Medium.Density rho[n] = Medium.density(state);
SI.InternalEnergy U[n] = m .* u;
Medium.SpecificInternalEnergy u[n] = Medium.specificInternalEnergy(state);
Medium.Temperature T[n] = Medium.temperature(state);
Medium.DynamicViscosity mu[n] = Medium.dynamicViscosity(state);
SI.Velocity v[n](nominal=0.2) = 0.5 * (mflow[1:n] + mflow[2:n+1]) ./ rho ./ A;
SI.Power Wflow[n];
SI.MomentumFlux Iflow[n] = v .* v .* rho * A;
// Variables faces: positiv in direction a -> b
Medium.MassFlowRate mflow[n+1](each start = mflow_start, each stateSelect=StateSelect.prefer, nominal=0.25) annotation(Dialog(tab="Initialization", showStartAttribute=true, enable=false));
Medium.EnthalpyFlowRate Hflow[n+1];
SI.Momentum I[n-1] = mflow[2:n] * dL;
SI.Force Fp[n-1];
SI.Force Ff[n-1];
SI.PressureDifference dpf[n-1](each start = (p_a_start - p_b_start)/(n-1), nominal=0.01e5) annotation(Dialog(tab="Initialization", showStartAttribute=true, enable=false));
equation
der(m) = mflow[1:n] - mflow[2:n+1]; // Mass balance
der(U) = Hflow[1:n] - Hflow[2:n+1] + Wflow; // Energy balance
der(I) = Iflow[1:n-1] - Iflow[2:n] + Fp - Ff; // Momentum balance, staggered
Hflow[1] = semiLinear(mflow[1], inStream(port_a.h_outflow), h[1]);
Hflow[2:n] = semiLinear(mflow[2:n], h[1:n-1], h[2:n]);
Hflow[n+1] = semiLinear(mflow[n+1], h[n], inStream(port_b.h_outflow));
Wflow[1] = v[1] * A .* ( (p[2] - p[1])/2 + dpf[1]/2);
Wflow[2:n-1] = v[2:n-1] * A .* ( (p[3:n]-p[1:n-2])/2 + (dpf[1:n-2]+dpf[2:n-1])/2);
Wflow[n] = v[n] * A .* ( (p[n] - p[n-1])/2 + dpf[n-1]/2);
Fp = A * (p[1:n-1] - p[2:n]);
Ff = A * dpf; // dpf = Ff ./ A;
if use_fixed_zeta then
dpf = 1/2 * zeta/(n-1) * (mflow[2:n]).^2 ./ ( 0.5*(rho[1:n-1] + rho[2:n]) * A * A);
else
dpf = homotopy(
actual = Modelica.Fluid.Pipes.BaseClasses.WallFriction.Detailed.pressureLoss_m_flow(
m_flow = mflow[2:n],
rho_a = rho[1:n-1],
rho_b = rho[2:n],
mu_a = mu[1:n-1],
mu_b = mu[2:n],
length = dL,
diameter = D,
roughness = R,
m_flow_small = 0.001),
simplified = dp_nominal/(n-1)/m_flow_nominal*mflow[2:n]);
end if;
// Boundary conditions
mflow[1] = port_a.m_flow;
mflow[n] = -port_b.m_flow;
p[1] = port_a.p;
p[n] = port_b.p;
port_a.h_outflow = h[1];
port_b.h_outflow = h[n];
initial equation
der(mflow[2:n]) = zeros(n-1);
der(p) = zeros(n);
der(h) = zeros(n);
annotation (Icon(coordinateSystem(preserveAspectRatio=false), graphics={Rectangle(
extent={{-100,60},{100,-60}},
fillColor={255,255,255},
fillPattern=FillPattern.HorizontalCylinder,
lineColor={0,0,0}),
Line(
points={{-100,60},{-100,-60}},
color={0,0,0},
thickness=0.5),
Line(
points={{-60,60},{-60,-60}},
color={0,0,0},
thickness=0.5),
Line(
points={{-20,60},{-20,-60}},
color={0,0,0},
thickness=0.5),
Line(
points={{20,60},{20,-60}},
color={0,0,0},
thickness=0.5),
Line(
points={{60,60},{60,-60}},
color={0,0,0},
thickness=0.5),
Line(
points={{100,60},{100,-60}},
color={0,0,0},
thickness=0.5),
Line(
points={{60,-80},{-60,-80}},
color={0,128,255},
visible=showDesignFlowDirection),
Polygon(
points={{20,-65},{60,-80},{20,-95},{20,-65}},
lineColor={0,128,255},
fillColor={0,128,255},
fillPattern=FillPattern.Solid,
visible=showDesignFlowDirection),
Text(
extent={{-150,100},{150,60}},
lineColor={0,0,255},
textString="%name"),
Text(
extent={{-40,22},{40,-18}},
lineColor={0,0,0},
textString="n = %n")}), Diagram(
coordinateSystem(preserveAspectRatio=false)));
end Pipe_FVM_staggered;
很长一段时间以来,我一直在努力解决这个问题,所以非常感谢您的任何评论或提示!!
如果您需要更多信息或测试结果,请告诉我!
这是
:
model Test_Waterhammer
extends Modelica.Icons.Example;
import SI = Modelica.SIunits;
import g = Modelica.Constants.g_n;
replaceable package Medium = Modelica.Media.Water.StandardWater;
Modelica.Fluid.Sources.Boundary_pT outlet(
redeclare package Medium = Medium,
nPorts=1,
p=2000000,
T=293.15)
annotation (Placement(transformation(extent={{90,-10},{70,10}})));
inner Modelica.Fluid.System system(
allowFlowReversal=true,
energyDynamics=Modelica.Fluid.Types.Dynamics.SteadyStateInitial,
massDynamics=Modelica.Fluid.Types.Dynamics.SteadyStateInitial,
momentumDynamics=Modelica.Fluid.Types.Dynamics.SteadyStateInitial,
m_flow_start=0.1,
m_flow_small=0.0001)
annotation (Placement(transformation(extent={{60,60},{80,80}})));
Modelica.Fluid.Sources.MassFlowSource_T inlet(
redeclare package Medium = Medium,
nPorts=1,
m_flow=0.1,
use_m_flow_in=true,
T=293.15)
annotation (Placement(transformation(extent={{-50,-10},{-30,10}})));
Modelica.Blocks.Sources.TimeTable timeTable(table=[0,0.1; 1,0.1; 1,0.25;
40,0.25; 40,0.35; 60,0.35])
annotation (Placement(transformation(extent={{-90,10},{-70,30}})));
Pipe_FVM_staggered pipe(
redeclare package Medium = Medium,
R=0.035*0.005,
mflow_start=0.1,
L=1000,
m_flow_nominal=0.1,
D=0.035,
zeta=2000,
n=11,
use_fixed_zeta=false,
T_start=293.15,
p_a_start=2010000,
p_b_start=2000000,
dp_nominal=10000)
annotation (Placement(transformation(extent={{10,-10},{30,10}})));
Modelica.Fluid.Pipes.DynamicPipe pipeMSL(
redeclare package Medium = Medium,
allowFlowReversal=true,
length=1000,
roughness=0.035*0.005,
m_flow_start=0.1,
energyDynamics=Modelica.Fluid.Types.Dynamics.SteadyStateInitial,
massDynamics=Modelica.Fluid.Types.Dynamics.SteadyStateInitial,
momentumDynamics=Modelica.Fluid.Types.Dynamics.SteadyStateInitial,
diameter=0.035,
modelStructure=Modelica.Fluid.Types.ModelStructure.av_vb,
redeclare model FlowModel =
Modelica.Fluid.Pipes.BaseClasses.FlowModels.DetailedPipeFlow (
useUpstreamScheme=false, use_Ib_flows=true),
p_a_start=2010000,
p_b_start=2000000,
T_start=293.15,
nNodes=11)
annotation (Placement(transformation(extent={{10,-50},{30,-30}})));
Modelica.Fluid.Sources.MassFlowSource_T inlet1(
redeclare package Medium = Medium,
nPorts=1,
m_flow=0.1,
use_m_flow_in=true,
T=293.15)
annotation (Placement(transformation(extent={{-48,-50},{-28,-30}})));
Modelica.Fluid.Sources.Boundary_pT outlet1(
redeclare package Medium = Medium,
nPorts=1,
p=2000000,
T=293.15)
annotation (Placement(transformation(extent={{90,-50},{70,-30}})));
equation
connect(inlet.ports[1], pipe.port_a)
annotation (Line(points={{-30,0},{-10,0},{10,0}}, color={0,127,255}));
connect(pipe.port_b, outlet.ports[1])
annotation (Line(points={{30,0},{50,0},{70,0}}, color={0,127,255}));
connect(inlet1.ports[1], pipeMSL.port_a)
annotation (Line(points={{-28,-40},{-10,-40},{10,-40}}, color={0,127,255}));
connect(pipeMSL.port_b, outlet1.ports[1])
annotation (Line(points={{30,-40},{50,-40},{70,-40}}, color={0,127,255}));
connect(timeTable.y, inlet.m_flow_in)
annotation (Line(points={{-69,20},{-60,20},{-60,8},{-50,8}}, color={0,0,127}));
connect(inlet1.m_flow_in, inlet.m_flow_in)
annotation (Line(points={{-48,-32},{-60,-32},{-60,8},{-50,8}}, color={0,0,127}));
annotation (Icon(coordinateSystem(preserveAspectRatio=false)), Diagram(
coordinateSystem(preserveAspectRatio=false)),
experiment(
StopTime=15,
__Dymola_NumberOfIntervals=6000,
Tolerance=1e-005,
__Dymola_Algorithm="Dassl"));
end Test_Waterhammer;
用301电池测试
:
解决方案:根据scottG的建议进行修改
model FVM_staggered_Ncells
// Import
import SI = Modelica.SIunits;
import Modelica.Constants.pi;
// Medium
replaceable package Medium = Modelica.Media.Interfaces.PartialMedium "Medium in the component"
annotation (choicesAllMatching = true);
// Interfaces, Ports
Modelica.Fluid.Interfaces.FluidPort_a port_a(redeclare package Medium = Medium) annotation (Placement(transformation(extent={{-110,-10},{-90,10}})));
Modelica.Fluid.Interfaces.FluidPort_b port_b(redeclare package Medium = Medium) annotation (Placement(transformation(extent={{90,-10},{110,10}})));
// Parameters
parameter Integer n(min=2) = 3 "Number of cells"; // No effect yet, only for icon
parameter SI.Length L = 1 "Length";
parameter SI.Diameter D = 0.010 "Diameter";
parameter SI.Height R = 2.5e-5 "Roughness";
parameter Boolean use_fixed_zeta = false "Use fixed zeta value instead of Moody chart";
parameter SI.CoefficientOfFriction zeta = 1;
// Initialization
parameter Medium.Temperature T_start = 293.15 "Start temperature" annotation(Dialog(tab="Initialization"));
parameter Medium.MassFlowRate mflow_start = 1 "Start mass flow rate in design direction" annotation(Dialog(tab="Initialization"));
parameter Medium.AbsolutePressure p_a_start = 2e5 "Start pressure p[1] at design inflow" annotation(Dialog(tab="Initialization"));
parameter Medium.AbsolutePressure p_b_start = 1e5 "Start pressure for p[n+1] at design outflow" annotation(Dialog(tab="Initialization"));
parameter Medium.AbsolutePressure p_start[:] = linspace(p_a_start, p_b_start, n) annotation(Dialog(tab="Initialization"));
// parameter Medium.SpecificEnthalpy h_start[:] = Medium.specificEnthalpy_pTX(p_start, T_start, Medium.X_default);
parameter Medium.SpecificEnthalpy h_start = Medium.specificEnthalpy_pTX((p_a_start + p_b_start)/2, T_start, Medium.X_default) annotation(Dialog(tab="Initialization"));
parameter SI.AbsolutePressure dp_nominal = 1e5;
parameter SI.MassFlowRate m_flow_nominal = 1;
// Variables general
SI.Length dL = L/n;
SI.Length dLs[n-1] = cat(1,{1.5*dL}, fill(dL,n-3), {1.5*dL});
SI.Area A = D^2*pi/4;
SI.Volume V = A * dL;
// Variables cell centers: positiv in direction a -> b
Medium.AbsolutePressure p[n](start = p_start, each stateSelect=StateSelect.prefer) annotation(Dialog(tab="Initialization", showStartAttribute=true, enable=false));
Medium.SpecificEnthalpy h[n](each start = h_start, each stateSelect=StateSelect.prefer) annotation(Dialog(tab="Initialization", showStartAttribute=true, enable=false));
Medium.ThermodynamicState state[n] = Medium.setState_phX(p,h);
SI.Mass m[n] = rho .* V;
Medium.Density rho[n] = Medium.density(state);
SI.InternalEnergy U[n] = m .* u;
Medium.SpecificInternalEnergy u[n] = Medium.specificInternalEnergy(state);
Medium.Temperature T[n] = Medium.temperature(state);
Medium.DynamicViscosity mu[n] = Medium.dynamicViscosity(state);
SI.Velocity v[n] = 0.5 * (mflow[1:n] + mflow[2:n+1]) ./ rho ./ A;
SI.Power Wflow[n];
SI.MomentumFlux Iflow[n] = v .* v .* rho * A;
// Variables faces: positiv in direction a -> b
Medium.MassFlowRate mflow[n+1](each start = mflow_start, each stateSelect=StateSelect.prefer) annotation(Dialog(tab="Initialization", showStartAttribute=true, enable=false));
Medium.EnthalpyFlowRate Hflow[n+1];
SI.Momentum I[n-1] = mflow[2:n] .* dLs;
SI.Force Fp[n-1];
SI.Force Ff[n-1];
SI.PressureDifference dpf[n-1](each start = (p_a_start - p_b_start)/(n-1)) annotation(Dialog(tab="Initialization", showStartAttribute=true, enable=false));
equation
der(m) = mflow[1:n] - mflow[2:n+1]; // Mass balance
der(U) = Hflow[1:n] - Hflow[2:n+1] + Wflow; // Energy balance
der(I) = Iflow[1:n-1] - Iflow[2:n] + Fp - Ff; // Momentum balance, staggered
Hflow[1] = semiLinear(mflow[1], inStream(port_a.h_outflow), h[1]);
Hflow[2:n] = semiLinear(mflow[2:n], h[1:n-1], h[2:n]);
Hflow[n+1] = semiLinear(mflow[n+1], h[n], inStream(port_b.h_outflow));
Wflow[1] = v[1] * A .* ( (p[2] - p[1])/2 + dpf[1]/2);
Wflow[2:n-1] = v[2:n-1] * A .* ( (p[3:n]-p[1:n-2])/2 + (dpf[1:n-2]+dpf[2:n-1])/2);
Wflow[n] = v[n] * A .* ( (p[n] - p[n-1])/2 + dpf[n-1]/2);
Fp = A * (p[1:n-1] - p[2:n]);
Ff = A * dpf;
if use_fixed_zeta then
dpf = 0.5 * zeta/(n-1) * abs(mflow[2:n]) .* mflow[2:n] ./ ( 0.5*(rho[1:n-1] + rho[2:n]) * A * A);
else
dpf = homotopy(
actual = Modelica.Fluid.Pipes.BaseClasses.WallFriction.Detailed.pressureLoss_m_flow(
m_flow = mflow[2:n],
rho_a = 0.5*(rho[1:n-1] + rho[2:n]),
rho_b = 0.5*(rho[1:n-1] + rho[2:n]),
mu_a = 0.5*(mu[1:n-1] + mu[2:n]),
mu_b = 0.5*(mu[1:n-1] + mu[2:n]),
length = dLs,
diameter = D,
roughness = R,
m_flow_small = 0.001),
simplified = dp_nominal/(n-1)/m_flow_nominal*mflow[2:n]);
end if;
// Boundary conditions
mflow[1] = port_a.m_flow;
mflow[n+1] = -port_b.m_flow;
p[1] = port_a.p;
p[n] = port_b.p;
port_a.h_outflow = h[1];
port_b.h_outflow = h[n];
initial equation
der(mflow[2:n]) = zeros(n-1);
der(p) = zeros(n);
der(h) = zeros(n);
annotation (Icon(coordinateSystem(preserveAspectRatio=false), graphics={Rectangle(
extent={{-100,60},{100,-60}},
fillColor={255,255,255},
fillPattern=FillPattern.HorizontalCylinder,
lineColor={0,0,0}),
Line(
points={{-100,60},{-100,-60}},
color={0,0,0},
thickness=0.5),
Line(
points={{-60,60},{-60,-60}},
color={0,0,0},
thickness=0.5),
Line(
points={{-20,60},{-20,-60}},
color={0,0,0},
thickness=0.5),
Line(
points={{20,60},{20,-60}},
color={0,0,0},
thickness=0.5),
Line(
points={{60,60},{60,-60}},
color={0,0,0},
thickness=0.5),
Line(
points={{100,60},{100,-60}},
color={0,0,0},
thickness=0.5),
Line(
points={{60,-80},{-60,-80}},
color={0,128,255},
visible=showDesignFlowDirection),
Polygon(
points={{20,-65},{60,-80},{20,-95},{20,-65}},
lineColor={0,128,255},
fillColor={0,128,255},
fillPattern=FillPattern.Solid,
visible=showDesignFlowDirection),
Text(
extent={{-150,100},{150,60}},
lineColor={0,0,255},
textString="%name"),
Text(
extent={{-40,22},{40,-18}},
lineColor={0,0,0},
textString="n = %n")}),
Diagram(coordinateSystem(preserveAspectRatio=false)));
end FVM_staggered_Ncells;
