Pipe

model Pipe
    import HydraulicsByFluidon.Media.Base.FluidInterface;
    import SI = Modelica.SIunits;
    import HydraulicsByFluidon.Tools.Pipes;
    import HydraulicsByFluidon.Types.TFluidPort;

    parameter SI.Length diameter = 0.032 "inner diameter"
        annotation (Dialog(group = "Geometry"));
    parameter SI.Length length = 1 "pipe length"
        annotation (Dialog(group = "Geometry"));
    parameter SI.Length deltaH = 0 "height difference between the ports B and A"
        annotation (Dialog(group = "Geometry"));
    parameter Modelica.SIunits.DimensionlessRatio relRoughness = 1e-6 "inner surface relative roughness"
        annotation (Dialog(group = "Geometry"));
    parameter Boolean flowInitialisation = true "Initialization of flow rate"
        annotation (Dialog(group = "Initialization"));
    parameter SI.Frequency fMax = 10 "Highest frequency to be covered"
        annotation (Dialog(group = "Initialization"));
    parameter SI.Velocity cMin = 1000 "Lowest speed of sound to be expected"
        annotation (Dialog(group = "Initialization"));
    parameter Boolean forwardFluidProperties = true "Forward fluid properties between ports"
        annotation (Dialog(tab = "Fluid Properties"));
    parameter TFluidPort propertyPort = TFluidPort.portA "Property determining port if forwardFluidProperties = false"
        annotation (Dialog(tab = "Fluid Properties"));
    HydraulicsByFluidon.Components.Lines.PipeWithoutCapacity pipeN[nDynPipes](each forwardFluidProperties = true, each diameter = diameter, each length = length / (nDynPipes + 2), each deltaH = deltaH / (nDynPipes + 2), each relRoughness = relRoughness);
    HydraulicsByFluidon.Components.Lines.PipeWithoutCapacity pipeFirst(forwardFluidProperties = forwardFluidProperties or not forwardFluidProperties and propertyPort == TFluidPort.portA, diameter = diameter, length = length / (nDynPipes + 2), deltaH = deltaH / (nDynPipes + 2), relRoughness = relRoughness) annotation (Placement(
        visible = true,
        transformation(
            origin = {0, 0},
            extent = {
                {-19, -8}, 
                {21, 12}},
            rotation = 0)));
    HydraulicsByFluidon.Components.Lines.PipeWithoutCapacity pipeLast(forwardFluidProperties = forwardFluidProperties or not forwardFluidProperties and propertyPort == TFluidPort.portB, diameter = diameter, length = length / (nDynPipes + 2), deltaH = deltaH / (nDynPipes + 2), relRoughness = relRoughness);
    HydraulicsByFluidon.Components.Base.Capacity volumeN[nDynPipes + 1](each capacity = pipeVolume / (nDynPipes + 1));
    HydraulicsByFluidon.Interfaces.FluidPort fluidPortA(p(nominal = 100000)) "Hydraulic port A"
        annotation (Placement(
            visible = true,
            transformation(
                extent = {
                    {-10, -90}, 
                    {10, -110}},
                rotation = 0),
            iconTransformation(
                extent = {
                    {-10, -90}, 
                    {10, -110}},
                rotation = 0)));
    HydraulicsByFluidon.Interfaces.FluidPort fluidPortB(p(nominal = 100000)) "Hydraulic port B"
        annotation (Placement(
            visible = true,
            transformation(
                extent = {
                    {-10, 110}, 
                    {10, 90}},
                rotation = 0),
            iconTransformation(
                extent = {
                    {-10, 110}, 
                    {10, 90}},
                rotation = 0)));
protected
    outer HydraulicsByFluidon.Media.Environment environment;
    parameter SI.Volume pipeVolume = Modelica.Constants.pi * diameter ^ 2 * length * (0.25);
    parameter Integer nDynPipes = if 0 < nPipesTheo - 2 then integer(nPipesTheo) - 2 else 0;
    parameter Real iReal[nDynPipes + 1](each fixed = false);
    parameter Real nPipesTheo = ceil(10 * fMax * length / cMin);
    parameter SI.MassFlowRate mFlowInit(fixed = false);
    parameter SI.AbsolutePressure dpInit(fixed = false);
    parameter Integer i(fixed = false, start = 0);
    parameter Real dpOffset(fixed = false);
    parameter Real rho(fixed = false);
    parameter Real eta(fixed = false);

    function calcMFlowInit
        input Boolean flowInitialisation;
        input Real relRoughness;
        input Real dpInit;
        input Real rho;
        input Real eta;
        input Real diameter;
        input Real length;
        output Real mFlowInit;
    protected
        Real reInitNew;
        Real reInitOld = 2299;
        Real lambdaInit;
        Integer i = 0;
        Real area = Modelica.Constants.pi * diameter ^ 2 * (0.25);
    algorithm
        if flowInitialisation and 0 < abs(dpInit) then 
            reInitNew := abs(dpInit) * rho * diameter ^ 3 / (32 * length * eta ^ 2);
            if 2320 < reInitNew then 
                while i < 25 and 1e-9 < abs(reInitNew / reInitOld - 1) loop
                    i := i + 1;
                    lambdaInit := HydraulicsByFluidon.Tools.Pipes.calcLambda(relRoughness, reInitNew);
                    mFlowInit := sign(dpInit) * area * sqrt(2 * rho * diameter * abs(dpInit) / (length * lambdaInit));
                    reInitOld := reInitNew;
                    reInitNew := mFlowInit * diameter / (eta * area);
                end while;
            end if;
            mFlowInit := 0.25 * sign(dpInit) * reInitNew * eta * Modelica.Constants.pi * diameter;
        else 
            mFlowInit := 0;
        end if;
    end calcMFlowInit;
initial algorithm
    rho := 0.5 * (FluidInterface.calcRho(fluidPortA.fluidId, fluidPortA.p, fluidPortA.fluidTemperature) + FluidInterface.calcRho(fluidPortB.fluidId, fluidPortB.p, fluidPortB.fluidTemperature));
    eta := rho * (0.5) * (FluidInterface.calcNu(fluidPortA.fluidId, fluidPortA.p, fluidPortA.fluidTemperature) + FluidInterface.calcNu(fluidPortB.fluidId, fluidPortB.p, fluidPortB.fluidTemperature));
    dpOffset := deltaH * environment.g * rho;
    dpInit := fluidPortA.p - fluidPortB.p - dpOffset;
    mFlowInit := calcMFlowInit(flowInitialisation, relRoughness, dpInit, rho, eta, diameter, length);
initial equation
    for i in 1:nDynPipes + 1 loop
        iReal[i] = i;
        volumeN[i].fluidPort.p = fluidPortA.p + (fluidPortB.p - fluidPortA.p) * iReal[i] / (nDynPipes + 2);
    end for;
    for i in 1:nDynPipes loop
        pipeN[i].fluidPortA.mFlow = mFlowInit;
    end for;
    assert(0 <= relRoughness, "Relative roughness muss be greater than 0!");
    assert(deltaH <= length, "Decrease deltaH or increase length!");
    assert(0.1 < cMin, "Increase cMin!");
    assert(0.01 < fMax, "Increase fMax!");
    assert(1e-15 < diameter, "Increase diameter!");
    assert(1e-15 < length, "Increase pipe length!");
    assert(10 * diameter < length, "Length has to be larger than 10 * diameter!");
    assert(nDynPipes < 18, "Too many pipe elements! Reduce fMax or increase cMin.");
    Modelica.Utilities.Streams.print("Number of internal pipes = " + String(nDynPipes + 2));
    pipeFirst.fluidPortA.mFlow = mFlowInit;
    pipeLast.fluidPortA.mFlow = mFlowInit;
equation
    connect(fluidPortA,pipeFirst.fluidPortA);
    connect(pipeFirst.fluidPortB,volumeN[1].fluidPort);
    connect(volumeN[nDynPipes + 1].fluidPort,pipeLast.fluidPortA);
    if 1 <= nDynPipes then 
        for i in 1:nDynPipes loop
            connect(volumeN[i].fluidPort,pipeN[i].fluidPortA);
            connect(pipeN[i].fluidPortB,volumeN[i + 1].fluidPort);
        end for;
    end if;
    connect(pipeLast.fluidPortB,fluidPortB) annotation (Line(
        points = {
            {0, 100}, 
            {0, 40}},
        color = {0, 93, 152}));

    annotation (
        defaultComponentName = "pipe",
        Icon(
            coordinateSystem(initialScale = 0.1),
            graphics = {
                Text(
                    origin = {200, -40},
                    lineColor = {0, 0, 255},
                    extent = {
                        {-160, 120}, 
                        {175, 80}},
                    textString = "%name",
                    horizontalAlignment = TextAlignment.Left), 
                Line(
                    origin = {-78, 104},
                    points = {
                        {0, 0}}), 
                Polygon(
                    origin = {21, 68},
                    lineColor = {0, 93, 152},
                    fillColor = {255, 255, 255},
                    fillPattern = FillPattern.VerticalCylinder,
                    points = {
                        {-39, 20}, 
                        {-45, 18}, 
                        {-50, 14.55}, 
                        {-51, 12}, 
                        {-51, -148}, 
                        {9, -148}, 
                        {9, 12}, 
                        {8, 14.55}, 
                        {3, 18}, 
                        {-3, 20}, 
                        {-15, 22}, 
                        {-27, 22}, 
                        {-39, 20}}), 
                Ellipse(
                    origin = {0, -87},
                    fillColor = {200, 200, 200},
                    fillPattern = FillPattern.VerticalCylinder,
                    extent = {
                        {-30, 17}, 
                        {30, -3}},
                    endAngle = 360), 
                Text(
                    origin = {0, -40},
                    lineColor = {0, 0, 255},
                    fillColor = {0, 0, 255},
                    extent = {
                        {-78, -20}, 
                        {78, 20}},
                    textString = "A"), 
                Text(
                    origin = {0, 60},
                    lineColor = {0, 0, 255},
                    fillColor = {0, 0, 255},
                    extent = {
                        {-78, -20}, 
                        {78, 20}},
                    textString = "B")}),
        Documentation(info = "<html>\n            <p>\n                Pipe model with resistive, capacitive and inductive properties. Depending on <var>pipe length</var>,\n                the <var>Highest frequency to be covered</var> and the <var>Lowest speed of sound to be expected</var> \n                during the simulation, the pipe model is subdivided into several sub-pipes\n                (inductive & resistive properties) and volumes (capacitive properties). Estimation of the required\n                number of pipes is based on the following equation:\n            </p>\n            <p>\n                <center><img align=\"middle\" src=\"modelica://HydraulicsByFluidon/Resources/Images/Components/Lines/EqLinesNoOfPipes.gif\"></center>\n            </p>\n            <p>\n                Irrespective of the parameterisation, at least two sub-pipes and one intermediate volume\n                will be created. The total volume of the pipe is distributed equally over the individual volumes\n            </p>\n            <p>\n                The pressure drop calculation is based on the assumption \n                of quasi-steady friction. For turbulent flow through hydraulically \n                rough pipes, the relative surface roughness (relRoughness) is taken into account. \n                Geodetic pressure offset due to height differences between the pipe\n                ports is considered in the simulation. <var>height difference between the ports B and A</var> is measured as the \n                vertical distance between port B and port A of the pipe.\n            </p>\n            <p>\n                The calculated inductance depends on the flow regime (laminar or turbulent). The laminar inductance\n                exceeds the turbulent one by a factor of 4/3.\n            </p>\n            <p>\n                If the parameter <var>Initialization of flow rate</var> is activated, the mass flow rate will be initialized based on\n                the potential difference between the two ports. If the check mark is not set, the initial flow rate will\n                be zero even if a potential difference is present.\n            </p></html>"));
end Pipe;