modelica - Dynamic pipe model in MSL, Finite Volume Method -

i'm trying use modelica modeling of system composed of elastic pipes. now, i'm trying implement own dynamic pipe model (rigid, not yet elastic) using same approach (finite volume, staggered) in modelica.fluid library, of course not including options.

this model should more simple understand, it's flat model, not extending other classes. important because therefore colleagues can understand model without modelica knowhow , can convince them modelica adequate tool our purposes!

as test case use mass flow source step signal (waterhammer). model gives not same results modelica.fluid component. appreciate, if can me, understand what's happening!

the test system looks this:

the results 11 cells this:

as can see, pressure peak higher msl component , frequency/period not same. when chose more cells error gets smaller.

i'm quite sure i'm using same equations. cause of numerical reasons (i tryied using nominal values)? included own "fixed zeta" flow model modelica.fluid component can compare in case of fixed pressure loss coefficient zeta.

the code of pipe model quite short , nice if work this:

model pipe_fvm_staggered    // import   import si = modelica.siunits;   import modelica.constants.pi;    // medium   replaceable package medium = modelica.media.interfaces.partialmedium "medium in 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, 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] @ design inflow" annotation(dialog(tab="initialization"));   parameter medium.absolutepressure p_b_start = 1e5 "start pressure p[n+1] @ 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 = * dl;    // variables cell centers: positiv in direction -> 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 -> 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] * .* ( (p[2] - p[1])/2 + dpf[1]/2);   wflow[2:n-1] = v[2:n-1] * .* ( (p[3:n]-p[1:n-2])/2 + (dpf[1:n-2]+dpf[2:n-1])/2);   wflow[n] = v[n] * .* ( (p[n] - p[n-1])/2 + dpf[n-1]/2);    fp = * (p[1:n-1] - p[2:n]);   ff = * dpf; // dpf = ff ./ a;    if use_fixed_zeta     dpf = 1/2 * zeta/(n-1) * (mflow[2:n]).^2 ./ ( 0.5*(rho[1:n-1] + rho[2:n]) * * 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; 

i'm struggling problem since quite long time, comments or hints appreciated!! if need more information or test results, please tell me!

this code test example:

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; 

i've run test 301 cells:

zoom peak 1 , 2:

solution: modifications suggested scottg

model fvm_staggered_ncells    // import   import si = modelica.siunits;   import modelica.constants.pi;    // medium   replaceable package medium = modelica.media.interfaces.partialmedium "medium in 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, 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] @ design inflow" annotation(dialog(tab="initialization"));   parameter medium.absolutepressure p_b_start = 1e5 "start pressure p[n+1] @ 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 = d^2*pi/4;   si.volume v = * dl;    // variables cell centers: positiv in direction -> 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 -> 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] * .* ( (p[2] - p[1])/2 + dpf[1]/2);   wflow[2:n-1] = v[2:n-1] * .* ( (p[3:n]-p[1:n-2])/2 + (dpf[1:n-2]+dpf[2:n-1])/2);   wflow[n] = v[n] * .* ( (p[n] - p[n-1])/2 + dpf[n-1]/2);    fp = * (p[1:n-1] - p[2:n]);   ff = * dpf;    if use_fixed_zeta     dpf = 0.5 * zeta/(n-1) *  abs(mflow[2:n]) .* mflow[2:n] ./ ( 0.5*(rho[1:n-1] + rho[2:n]) * * 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; 

right result:

alrighty.. after digging figured out. below have shown "as received" code , edit below it. fixes all.

background, know there model structure very important. 1 modeled av_vb.

1. correct length of flow model

the variable dl (length of flow segments) different first , last volume of av_vb model structure. correction important case run.

add following modification:

// define variable si.length dls[n-1]; si.momentum i[n-1] = mflow[2:n] .* dls; // changed *dl .*dls  // add equation section dls[1] = dl + 0.5*dl; dls[2:n-2] = fill(dl,n-3); dls[n-1] =  dl + 0.5*dl; 

2. change dpf mflow calculation

i ran simple case constant flow calculation , checked results , found different first correction. appears dpf = f(mflow) calculation used when under specified settings "one-to-one" comparison use mflow=f(dpf). because have chosen momentumdynamics=steadystateinitial makes from_dp = true in partialgenericpipeflow. if change results same constant flow example (differences between 2 show easier since not covered time dependent dynamics of changing flow rates).

also, average densities being used differed msl pipe think. didn't impact calculation example feel free double check conclusion.

  if use_fixed_zeta     dpf = 1/2*zeta/(n - 1)*(mflow[2:n]) .^ 2 ./ (0.5*(rho[1:n - 1] + rho[2:n])*       a*a);   else  // original     //      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 = dls, //notice changed dl dls     //          diameter = d,     //          roughness = r,     //          m_flow_small = 0.001),     //        simplified = dp_nominal/(n-1)/m_flow_nominal*mflow[2:n]);  // correct model "one-to-one" comparison chosen conditions. averaged rho , mu used since useupstreamscheme = false.     mflow[2:n] = homotopy(actual=       modelica.fluid.pipes.baseclasses.wallfriction.detailed.massflowrate_dp(       dpf,       0.5*(rho[1:n - 1] + rho[2:n]),       0.5*(rho[1:n - 1] + rho[2:n]),       0.5*(mu[1:n - 1] + mu[2:n]),       0.5*(mu[1:n - 1] + mu[2:n]),       dls,       d,       a,       r,       1e-5,       4000), simplified=m_flow_nominal/dp_nominal .* dpf);   end if; 

3. correct port_b.m_flow reference

this edit did not impact results of calculation in others.

// original   mflow[n] = -port_b.m_flow; // fixed reference proper flow variable   mflow[n+1] = -port_b.m_flow; 

below same plot generated. plots overlap.