Ticket #134: pipe.mso

using "streams";

Model props

ATTRIBUTES
        Pallete         = false;
        Brief           = "System properties for the pipe model";
        Info                    =
"== Contents ==
* Mass Density
* Dynamic Viscosity
* Molar Volume
* Reynolds Number
* Velocity
";

PARAMETERS

outer   N as Integer    (Brief = "Number of Profile Intervals", Default = 1, Lower = 1, Upper = 100);

VARIABLES

        rho(N+1)                as dens_mass    (Brief = "Mass Density Profile" , Symbol = "\rho");
        mu(N+1)                 as viscosity                    (Brief = "Viscosity Profile" , Symbol = "\mu");
        Re(N+1)                 as Real                         (Brief = "Reynolds Number Profile");
        Vel(N+1)                as velocity                     (Brief = "Velocity Profile");
        vm(N+1)                 as vol_mol                      (Brief = "Mixture Molar Volume Profile");
        
end

Model pipe

ATTRIBUTES
        Pallete         = true;
        Icon            = "icon/pipe";          
        Brief           = "pipe";
        Info            =
"This distributed model describes the pressure drop of a liquid stream flowing in a pipe.

==Assumptions==
*The flow along the pipe is adiabatic;
*The pipe is always full of liquid;
*Cross sectional area is constant;
*Newtonian fluid;
*Steady-state;
";

PARAMETERS

        outer   NComp   as Integer      (Brief = "Number of chemical components", Lower = 1); 
        outer  PP               as Plugin       (Brief = "External Physical Properties",Type="PP");

        N                               as Integer                      (Brief = "Number of Profile Intervals", Default = 1, Lower = 1, Upper = 100);
        pi                              as Real                                 (Brief="pi number",Default=3.141592, Symbol = "\pi");
        g                                       as acceleration         (Brief="Acceleration of gravity");
        Lpipe                   as length                       (Brief="Pipe Length", Symbol = "L_{pipe}");
        Hrise                   as length                       (Brief="Pipe Rise", Symbol = "H_{rise}");
        Dpipe                   as length                       (Brief="Pipe Inner Diameter", Symbol = "D_{pipe}");
        Apipe                   as area                         (Brief="Pipe Area", Symbol = "A_{pipe}");
        Roughness as length                     (Brief="Pipe Roughness", Symbol = "\varepsilon");

        FlowRegime      as Switcher     (Brief="Pipe flow regime ",Valid=["laminar","turbulent"],Default="laminar");

SET

        g                       = 1*'ga';
        Apipe   = 0.25*pi*Dpipe^2; 

VARIABLES
        
        in              Inlet                           as stream               (Brief = "Inlet Stream" ,PosX=0, PosY=0.5225, Symbol = "^{in}");
        out     Outlet                  as streamPH     (Brief = "Outlet Stream",PosX=1, PosY=0.5225, Symbol = "^{out}");
                        Properties      as props                        (Brief = "Pipe Properties", Symbol = " ");



        Pdrop                           as pressure     (Brief = "Total Pressure Drop", DisplayUnit = 'kPa',Lower = 0, Symbol = "\Delta P_{drop}");
        dPfric(N+1)     as pressure     (Brief = "Friction Pressure Drop", DisplayUnit = 'kPa',Lower = 0, Symbol = "\Delta P_{fric}");
        dPelv(N+1)      as pressure     (Brief = "Elevation Pressure Drop", DisplayUnit = 'kPa',Lower = 0 , Symbol = "\Delta P_{elev}");
        Pincr(N+1)              as pressure     (Brief = "Pressure Profile", DisplayUnit = 'kPa' , Symbol = "P_{incr}");
        Lincr(N+1)              as length               (Brief = "Length Points", Symbol = "L_{incr}");
        f(N+1)                  as fricfactor   (Brief = "Friction Factor");
        
EQUATIONS

"Inlet Boudary for Pressure Profile"
        Pincr(1) =  Inlet.P;

"Outlet Boundary for Pressure Profile"
        Pincr(N+1) = Outlet.P;

"Total Pressure Drop"
        Pdrop = dPfric(N+1) + dPelv(N+1);
        
"Pipe Initial Length"
        Lincr(1) = 0*'m';

"Outlet Composition"
        Outlet.z = Inlet.z;
        
"Outlet Temperature"
        Outlet.T = Inlet.T;
        
"Molar Balance"
        Outlet.F = Inlet.F;

"Velocity"
        Properties.Vel(1:N+1) = Inlet.F/Apipe*Properties.vm(1:N+1);

"Reynolds Number"
        Properties.Re(1:N+1)    = Properties.rho(1:N+1)*Properties.Vel(1:N+1)*Dpipe/Properties.mu(1:N+1);

"Incremental Friction Pressure Drop"
        dPfric(1:N+1)   = (2*f(1:N+1)*Lincr(1:N+1)*Properties.rho(1:N+1)*Properties.Vel(1:N+1)^2/Dpipe);

if Hrise > Lpipe
        
        then
"Incremental Elevation Pressure Drop Constraint"
        dPelv(1:N+1)    = Properties.rho(1:N+1)*g*Lincr(1:N+1)*1;
        
        else
"Incremental Elevation Pressure Drop"
        dPelv(1:N+1)    = Properties.rho(1:N+1)*g*Lincr(1:N+1)*(Hrise/abs(Lpipe));
        
end

for i in [1:N+1]

"Density"
        Properties.rho(i) = PP.LiquidDensity(Inlet.T,Pincr(i),Inlet.z);
        
"Viscosiyty"
        Properties.mu(i) = PP.LiquidViscosity(Inlet.T,Pincr(i),Inlet.z);
        
"Molar Volume"
        Properties.vm(i) = PP.LiquidVolume(Inlet.T,Pincr(i),Inlet.z);

end

for i in [1:N]
        
"Outlet Pressure"
        Pincr(i+1) = Pincr(1) - (dPfric(i+1) + dPelv(i+1));

"Incremental Length"
        Lincr(i+1) = i*abs(Lpipe)/N;

end

for i in [1:N+1]
        
switch FlowRegime
        
        case "laminar":
        
"Friction Factor for Pressure Drop - laminar Flow"
        f(i)*Properties.Re(i) = 16;
        
        when Properties.Re(i) > 2300 switchto "turbulent";

        case "turbulent":

"Friction Factor  for Pressure Drop - Turbulent Flow"
        1/sqrt(f(i))= -4*log(Roughness/Dpipe/3.7+1.255/Properties.Re(i)/sqrt(f(i)));

        when Properties.Re(i) <= 2300 switchto "laminar";
        
end

end

end



FlowSheet Pipe_simples

PARAMETERS
        PP      as Plugin(Brief="Physical Properties",
                Type="PP",
                Components = [ "water"],
                LiquidModel = "PR",
                VapourModel = "PR"
        );
        NComp   as Integer;

DEVICES
        Tube as pipe;
        Feed as simple_source;
        
SET
        NComp  = PP.NumberOfComponents;
        Tube.N  = 10;
        
CONNECTIONS
        Feed.Outlet to Tube.Inlet;

SPECIFY
        Feed.F = 10 * 'mol/s';
        Feed.P = 150 * 'kPa';
        Feed.T = 281.75 * 'K';
        Feed.MolarComposition = [1];

SET

        Tube.Lpipe                      = 30*'m';
        Tube.Hrise                      =       5*'m';
        
        Tube.Dpipe                      = 3*'in';
        Tube.Roughness  =4.572e-05*'m';
        
OPTIONS
        Dynamic = false;
        NLASolver(
                RelativeAccuracy = 1e-6
        );
end