Changeset 504 for branches


Ignore:
Timestamp:
Apr 16, 2008, 7:50:32 PM (14 years ago)
Author:
gerson bicca
Message:

updated pipe model (improved model to calculate the pressure drop for two phase flow)

Location:
branches/tests
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed
  • branches/tests/eml/pressure_changers/pipe.mso

    r343 r504  
    1 #*-------------------------------------------------------------------
    2 * EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
    3 *
    4 * This LIBRARY is free software; you can distribute it and/or modify
    5 * it under the therms of the ALSOC FREE LICENSE as available at
    6 * http://www.enq.ufrgs.br/alsoc.
    7 *
    8 * EMSO Copyright (C) 2004 - 2007 ALSOC, original code
    9 * from http://www.rps.eng.br Copyright (C) 2002-2004.
    10 * All rights reserved.
    11 *
    12 * EMSO is distributed under the therms of the ALSOC LICENSE as
    13 * available at http://www.enq.ufrgs.br/alsoc.
    14 *----------------------------------------------------------------------
    15 * Author: Gerson Balbueno Bicca
    16 * $Id: pipe.mso 326 2007-08-17 19:25:59Z bicca $
    17 *--------------------------------------------------------------------*#
     1
    182using "streams";
     3
     4Model props
     5
     6ATTRIBUTES
     7        Pallete         = false;
     8        Brief           = "System properties for the pipe model";
     9
     10PARAMETERS
     11
     12outer   N                               as Integer      (Brief = "Number of Profile Intervals", Default = 1, Lower = 1, Upper = 100);
     13outer   NComp   as Integer      (Brief = "Number of chemical components", Lower = 1);
     14       
     15VARIABLES
     16
     17        Fw                                                      as flow_mass            (Brief = "Mass Flow Profile" , Symbol = "F_w");
     18        rho(N+1)                                as dens_mass            (Brief = "Mass Density Profile" , Symbol = "\rho");
     19        mu(N+1)                                 as viscosity                    (Brief = "Viscosity Profile" , Symbol = "\mu");
     20        Re(N+1)                                         as Real                                 (Brief = "Reynolds Number Profile");
     21        Vel(N+1)                                        as velocity                     (Brief = "Velocity Profile");
     22        vm(N+1)                                 as vol_mol                      (Brief = "Mixture Molar Volume Profile");
     23        Lfrac(N+1)                              as fraction                     (Brief = "Liquid Fraction Profile");
     24        Vfrac(N+1)                              as fraction                     (Brief = "Vapour Fraction Profile");
     25        x(N+1,NComp)            as fraction                     (Brief = "Liquid Molar Fraction");
     26        y(N+1,NComp)            as fraction                     (Brief = "Vapour Molar Fraction");
     27        h(N+1)                                          as enth_mol             (Brief = "Stream Enthalpy profile");
     28       
     29end
    1930
    2031Model pipe
     
    2334        Pallete         = true;
    2435        Icon            = "icon/pipe";         
    25         Brief           = "Model of a simple pipe";
     36        Brief           = "pipe";
    2637        Info            =
    27         "testing a model of pressure drop.";
     38"This distributed model describes the pressure drop of a material stream flowing in a pipe.
     39
     40==Assumptions==
     41*Cross sectional area is constant;
     42*Newtonian fluid;
     43*Steady-state;
     44";
    2845
    2946PARAMETERS
    3047
    31         outer   NComp   as Integer(Brief = "Number of chemical components", Lower = 1);
    32         outer PP                as Plugin (Brief = "External Physical Properties",Type="PP");
    33 
    34         N                                       as Integer                      (Brief = "Number of Profile Intervals", Default = 1, Lower = 1, Upper = 300);
    35         pi                              as Real                                 (Brief="pi number",Default=3.141592);
     48        outer   NComp   as Integer      (Brief = "Number of chemical components", Lower = 1);
     49        outer  PP               as Plugin       (Brief = "External Physical Properties",Type="PP");
     50
     51        N                               as Integer                      (Brief = "Number of Profile Intervals", Default = 1, Lower = 1, Upper = 100);
     52        pi                              as Real                                 (Brief="pi number",Default=3.141592, Symbol = "\pi");
    3653        g                                       as acceleration         (Brief="Acceleration of gravity");
    37         Lpipe                   as length                       (Brief="Pipe Length");
    38         Hrise                           as length                       (Brief="Pipe Rise");
    39         Dpipe                   as length                       (Brief="Pipe Inner Diameter");
    40         Apipe                   as area                         (Brief="Pipe Area");
    41         Roughness               as length                       (Brief="Pipe Roughness");
    42 
    43         FlowRegime      as Switcher     (Brief="Pipe flow regime ",Valid=["laminar","turbulent"],Default="laminar");
     54        Lpipe                   as length                       (Brief="Pipe Length", Symbol = "L_{pipe}");
     55        Hrise                   as length                       (Brief="Pipe Rise", Symbol = "H_{rise}");
     56        Dpipe                   as length                       (Brief="Pipe Inner Diameter", Symbol = "D_{pipe}");
     57        Apipe                   as area                         (Brief="Pipe Area", Symbol = "A_{pipe}");
     58        Roughness as length                             (Brief="Pipe Roughness", Symbol = "\varepsilon");
     59        M(NComp)        as molweight    (Brief="Component Mol Weight");
     60        FlowRegime      as Switcher             (Brief="Pipe flow regime ",Valid=["laminar","turbulent"],Default="laminar");
     61        Thermal                 as Switcher             (Brief="Pipe Thermal Specification ",Valid=["Constant Temperature","Linear Temperature profile"],Default="Constant Temperature");
     62        Toutlet                         as temperature (Brief = "Outlet Temperature", Symbol = "T_{out}");
    4463
    4564SET
    4665
    47         g               = 1*'ga';
    48         Apipe = 0.25*pi*Dpipe^2;
    49 
    50 SUBMODELS
    51        
    52         in              Inlet           as stream               (Brief = "Inlet Stream" , Symbol = "^in");
    53         out     Outlet  as streamPH     (Brief = "Outlet Stream", Symbol = "^out");
    54 
     66        g                       = 1*'ga';
     67        Apipe   = 0.25*pi*Dpipe^2;
     68        M               = PP.MolecularWeight();
     69       
    5570VARIABLES
    56 
     71       
     72        in              Inlet                           as stream               (Brief = "Inlet Stream" ,PosX=0, PosY=0.5225, Symbol = "^{in}");
     73        out     Outlet                  as streamPH     (Brief = "Outlet Stream",PosX=1, PosY=0.5225, Symbol = "^{out}");
     74                        Properties     as props                 (Brief = "Pipe Properties", Symbol = " ");
     75        Tincr(N+1)                      as temperature (Brief = "Temperature Profile", Symbol = "T_{incr}");
     76       
    5777        Pdrop                   as pressure     (Brief = "Total Pressure Drop", DisplayUnit = 'kPa',Lower = 0, Symbol = "\Delta P_{drop}");
     78        dPfricTotal             as pressure     (Brief = "Total Friction Pressure Drop", DisplayUnit = 'kPa',Lower = 0, Symbol = "\Delta Ptotal_{fric}");
     79        dPelvTotal      as pressure     (Brief = "Total Elevation Pressure Drop", DisplayUnit = 'kPa',Lower = 0 , Symbol = "\Delta Ptotal_{elev}");     
     80        dPaccTotal      as pressure     (Brief = "Total Acceleration Pressure Drop", DisplayUnit = 'kPa',Lower = 0 , Symbol = "\Delta Ptotal_{acc}");   
    5881        dPfric(N+1)     as pressure     (Brief = "Friction Pressure Drop", DisplayUnit = 'kPa',Lower = 0, Symbol = "\Delta P_{fric}");
    5982        dPelv(N+1)      as pressure     (Brief = "Elevation Pressure Drop", DisplayUnit = 'kPa',Lower = 0 , Symbol = "\Delta P_{elev}");
    60         Pincr(N+1)              as pressure             (Brief = "Pressure Profile", DisplayUnit = 'kPa' , Symbol = "P_{incr}");
    61         Lincr(N+1)              as length               (Brief = "Length Points", Symbol = "L_{incr}");
    62         Vel(N+1)                as velocity             (Brief = "Velocity Profile");
    63         vm(N+1)                 as vol_mol              (Brief = "Mixture Molar Volume Profile");
     83        dPacc(N+1)      as pressure     (Brief = "Acceleration Pressure Drop", DisplayUnit = 'kPa',Lower = 0 , Symbol = "\Delta P_{acc}");
     84        Pincr(N+1)      as pressure     (Brief = "Pressure Profile", DisplayUnit = 'kPa' , Symbol = "P_{incr}");
     85        Lincr(N+1)      as length               (Brief = "Length Points", Symbol = "L_{incr}");
    6486        f(N+1)                  as fricfactor   (Brief = "Friction Factor");
    6587       
    66         rho(N+1)        as dens_mass    (Brief = "Mass Density Profile" , Symbol = "\rho");
    67         mu(N+1)         as viscosity    (Brief = "Viscosity Profile" , Symbol = "\mu");
    68         Re(N+1)         as Real                 (Brief = "Reynolds Number Profile");
    69 
    7088EQUATIONS
    7189
     90"Mass Flow"
     91        Properties.Fw = sum(M*Inlet.z)*Inlet.F;
     92
     93"Inlet Boudary for Temperature Profile"
     94        Tincr(1) =  Inlet.T;
     95
     96"Outlet Boundary for Temperature Profile"
     97        Tincr(N+1) = Outlet.T;
     98       
    7299"Inlet Boudary for Pressure Profile"
    73100        Pincr(1) =  Inlet.P;
     
    77104
    78105"Total Pressure Drop"
    79         Pdrop = dPfric(N+1) + dPelv(N+1);
    80        
     106        Pdrop = dPfricTotal+ dPelvTotal + dPaccTotal;
     107
     108"Total Pressure Drop for Friction"
     109        dPfricTotal = dPfric(N+1);
     110
     111"Total Pressure Drop for Elevation"
     112        dPelvTotal = dPelv(N+1);
     113
     114"Total Pressure Drop for Acceleration "
     115        dPaccTotal= sum(dPacc);
     116
    81117"Pipe Initial Length"
    82118        Lincr(1) = 0*'m';
     
    85121        Outlet.z = Inlet.z;
    86122       
    87 "Outlet Temperature"
    88         Outlet.T = Inlet.T;
    89        
    90123"Molar Balance"
    91124        Outlet.F = Inlet.F;
    92125
     126"Velocity"
     127        Properties.Vel(1:N+1) = Inlet.F/Apipe*Properties.vm(1:N+1);
     128
     129"Reynolds Number"
     130        Properties.Re(1:N+1)    = Properties.rho(1:N+1)*Properties.Vel(1:N+1)*Dpipe/Properties.mu(1:N+1);
     131
     132"Incremental Friction Pressure Drop"
     133        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);
     134
     135"Incremental Acceleration Pressure Drop at Pipe Entry"
     136        dPacc(1) = 0*'Pa';
     137
     138if Hrise > Lpipe
     139       
     140        then
     141"Incremental Elevation Pressure Drop Constraint"
     142        dPelv(1:N+1)    = Properties.rho(1:N+1)*g*Lincr(1:N+1);
     143       
     144        else
     145"Incremental Elevation Pressure Drop"
     146        dPelv(1:N+1)    = Properties.rho(1:N+1)*g*Lincr(1:N+1)*(Hrise/abs(Lpipe));
     147       
     148end
     149
    93150for i in [1:N+1]
    94151
    95 "Velocity"
    96         Vel(i) = Inlet.F/Apipe*vm(i);
    97 
    98 "Reynolds Number"
    99         Re(i)   = rho(i)*Vel(i)*Dpipe/mu(i);
     152"Flash Calculation"
     153        [Properties.Vfrac(i), Properties.x(i,:), Properties.y(i,:)] = PP.FlashPH(Pincr(i), Properties.h(i), Inlet.z);
     154
     155"Liquid Fraction"
     156        Properties.Lfrac(i) = 1-Properties.Vfrac(i);
     157
     158"Enthalpy"
     159        Properties.h(i) = (1-Properties.Vfrac(i))*PP.LiquidEnthalpy(Tincr(i),Pincr(i),Properties.x(i,:)) + Properties.Vfrac(i)*PP.VapourEnthalpy(Tincr(i),Pincr(i),Properties.y(i,:));
    100160
    101161"Density"
    102         rho(i) = PP.LiquidDensity(Inlet.T,Pincr(i),Inlet.z);
     162        Properties.rho(i) = (1-Properties.Vfrac(i))*PP.LiquidDensity(Tincr(i),Pincr(i),Properties.x(i,:))+Properties.Vfrac(i)*PP.VapourDensity(Tincr(i),Pincr(i),Properties.y(i,:));
    103163       
    104164"Viscosiyty"
    105         mu(i) = PP.LiquidViscosity(Inlet.T,Pincr(i),Inlet.z);
     165        Properties.mu(i) = (1-Properties.Vfrac(i))*PP.LiquidViscosity(Tincr(i),Pincr(i),Properties.x(i,:))+Properties.Vfrac(i)*PP.VapourViscosity(Tincr(i),Pincr(i),Properties.y(i,:));
    106166       
    107167"Molar Volume"
    108         vm(i) = PP.LiquidVolume(Inlet.T,Pincr(i),Inlet.z);
    109 
    110 "Incremental Friction Pressure Drop"
    111         dPfric(i)       = (2*f(i)*Lincr(i)*rho(i)*Vel(i)^2/Dpipe);
    112 
    113 "Incremental Elevation Pressure Drop"
    114         dPelv(i)        = rho(i)*g*Lincr(i)*(Hrise/Lpipe);
     168        Properties.vm(i) = (1-Properties.Vfrac(i))*PP.LiquidVolume(Tincr(i),Pincr(i),Properties.x(i,:))+Properties.Vfrac(i)*PP.VapourVolume(Tincr(i),Pincr(i),Properties.y(i,:));
    115169
    116170end
     
    119173       
    120174"Outlet Pressure"
    121         Pincr(i+1) = Pincr(1) - (dPfric(i+1) + dPelv(i+1));
    122 
     175        Pincr(i+1) = Pincr(1) - (dPfric(i+1) + dPelv(i+1) + dPacc(i+1));
     176
     177"Incremental Acceleration Pressure Drop"
     178        #dPacc(i+1) = (Properties.Fw/Apipe)^2*(1/Properties.rho(i+1)-1/Properties.rho(i));
     179        dPacc(i+1) = dPacc(1);
     180       
    123181"Incremental Length"
    124         Lincr(i+1) = i*Lpipe/N;
     182        Lincr(i+1) = i*abs(Lpipe)/N;
     183
     184end
     185
     186for  i in [1:N]
     187
     188switch Thermal
     189
     190        case "Constant Temperature":
     191                "Outlet Temperature"
     192                Tincr(i+1) = Inlet.T;
     193       
     194        case "Linear Temperature profile":
     195                "Outlet Temperature"
     196                Toutlet - Tincr(1) = Lpipe*((Tincr(i+1)-Tincr(i))/(Lincr(i+1)-Lincr(i)));
     197
     198end
    125199
    126200end
     
    133207       
    134208"Friction Factor for Pressure Drop - laminar Flow"
    135         f(i)*Re(i) = 16;
    136        
    137         when Re(i) > 2300 switchto "turbulent";
     209        f(i)*Properties.Re(i) = 16;
     210       
     211        when Properties.Re(i) > 2300 switchto "turbulent";
    138212
    139213        case "turbulent":
    140214
    141215"Friction Factor  for Pressure Drop - Turbulent Flow"
    142         1/sqrt(f(i))= -4*log(Roughness/Dpipe/3.7+1.255/Re(i)/sqrt(f(i)));
    143 
    144         when Re(i) <= 2300 switchto "laminar";
    145        
    146 end
    147 
    148 end
    149 
    150 end
     216        1/sqrt(f(i))= -4*log(Roughness/Dpipe/3.7+1.255/Properties.Re(i)/sqrt(f(i)));
     217
     218        when Properties.Re(i) <= 2300 switchto "laminar";
     219       
     220end
     221
     222end
     223
     224end
  • branches/tests/sample/pressure_changers/sample_pipe.mso

    r343 r504  
    1313* available at http://www.enq.ufrgs.br/alsoc.
    1414*
    15 * This sample file needs VRTherm DEMO (www.vrtech.com.br) to run.
     15* This sample file needs VRTherm (www.vrtech.com.br) to run.
    1616*
    1717*----------------------------------------------------------------------
     
    2828        PP      as Plugin(Brief="Physical Properties",
    2929                Type="PP",
    30                 Components = [ "water"],
     30                Components = [ "water" ],
    3131                LiquidModel     = "PR",
    3232                VapourModel     = "PR"
     
    8989        PP      as Plugin(Brief="Physical Properties",
    9090                Type="PP",
    91                 Components = [ "water"],
    92                 LiquidModel = "PR",
    93                 VapourModel = "PR"
     91                Components = [ "isobutane", "benzene"],
     92                LiquidModel     = "PR",
     93                VapourModel     = "PR"
    9494        );
     95
    9596        NComp   as Integer;
    9697
     
    9899        Tube as pipe;
    99100        Feed as simple_source;
     101        Storage as simple_sink;
    100102       
    101103SET
    102104        NComp  = PP.NumberOfComponents;
    103         Tube.N  = 100;
     105        Tube.N  = 10;
     106       
     107        #Tube.Thermal ="Linear Temperature profile";
     108        #Tube.Toutlet  = (90+273.15) * 'K';
     109       
     110        Tube.Thermal ="Constant Temperature";
    104111       
    105112CONNECTIONS
    106113        Feed.Outlet to Tube.Inlet;
     114        Tube.Outlet to Storage.Inlet;
    107115
    108116SPECIFY
    109         Feed.Outlet.F = 10 * 'mol/s';
    110         Feed.Outlet.P = 150 * 'kPa';
    111         Feed.Outlet.T = 281.75 * 'K';
    112         Feed.Outlet.z = [1];
     117        Feed.Outlet.F = 10 * 'kmol/h';
     118        Feed.Outlet.P = 3* 'atm';
     119        Feed.Outlet.T = (110+273.15) * 'K';
     120        Feed.Outlet.z = [0.5,0.5];
    113121
    114122SET
    115         Tube.Lpipe                      = 30*'m';
    116         Tube.Hrise                      =       5*'m';
     123        Tube.Lpipe                      = 5*'m';
     124        Tube.Hrise                      =       0*'m';
    117125       
    118126
     
    123131        OPTIONS
    124132        Dynamic = false;
    125         NLASolver(
    126                 RelativeAccuracy = 1e-6
    127         );
    128133end
Note: See TracChangeset for help on using the changeset viewer.