Changeset 531 for branches


Ignore:
Timestamp:
May 30, 2008, 8:44:00 PM (14 years ago)
Author:
gerson bicca
Message:

added Beggs and Brill method for two phase pipe (is not working for all cases)

Location:
branches/tests
Files:
2 deleted
1 edited

Legend:

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

    r504 r531  
    1515VARIABLES
    1616
    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");
     17        Fw(N+1)                                 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
     21        Vel(N+1)                                        as velocity     (Brief = "Velocity Profile");
     22        vm(N+1)                                 as vol_mol      (Brief = "Molar Volume Profile");
     23        frac(N+1)                                       as fraction             (Brief = "Molar Fraction Profile");
     24        Vsfrac(N+1)                             as fraction             (Brief = "No Slip Volume Fraction Profile",Lower=0);
     25        Holdup(N+1)                             as fraction             (Brief = "Holdup Profile",Lower=0);
     26        Mfrac(N+1,NComp)        as fraction             (Brief = "Phase Molar Fraction Profile");
     27        h(N+1)                                          as enth_mol     (Brief = "Molar Enthalpy profile");
     28       
    2829       
    2930end
     
    4950        outer  PP               as Plugin       (Brief = "External Physical Properties",Type="PP");
    5051
    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");
    53         g                                       as acceleration         (Brief="Acceleration of gravity");
    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}");
     52        N                                               as Integer                              (Brief = "Number of Profile Intervals", Default = 1, Lower = 1, Upper = 100);
     53        pi                                      as Real                                 (Brief="pi number",Default=3.141592, Symbol = "\pi");
     54        eps                             as Real                                 (Brief="Small Number",Default=1E-8);
     55        g                                               as acceleration                 (Brief="Acceleration of gravity");
     56        Lpipe                           as length                               (Brief="Pipe Length", Symbol = "L_{pipe}");
     57        Hrise                           as angle                                (Brief="Pipe Rise", Symbol = "H_{rise}");
     58        Dpipe                           as length                               (Brief="Pipe Inner Diameter", Symbol = "D_{pipe}");
     59        Apipe                           as area                                 (Brief="Pipe Area", Symbol = "A_{pipe}");
     60        Roughness               as length                               (Brief="Pipe Roughness", Symbol = "\varepsilon");
     61        M(NComp)                as molweight            (Brief="Component Mol Weight");
     62        FlowRegime              as Switcher                     (Brief="Pipe flow regime",Valid=["laminar","turbulent"],Default="laminar");
     63        Correlation             as Switcher                     (Brief="Holdup Correlation for Two Phase Flow", Valid=["Beggs-Brill","Lockhart-Martinelli"], Default="Beggs-Brill");
     64        Thermal                         as Switcher                     (Brief="Pipe Thermal Specification ",Valid=["Constant Temperature","Linear Temperature profile"],Default="Constant Temperature");
     65        Toutlet                                 as temperature          (Brief= "Outlet Temperature", Symbol = "T_{out}");
    6366
    6467SET
     
    7073VARIABLES
    7174       
    72         in              Inlet                           as stream               (Brief = "Inlet Stream" ,PosX=0, PosY=0.5225, Symbol = "^{in}");
     75        in              Inlet                   as stream               (Brief = "Inlet Stream" ,PosX=0, PosY=0.5225, Symbol = "^{in}");
    7376        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        
    77         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}");   
    81         dPfric(N+1)     as pressure     (Brief = "Friction Pressure Drop", DisplayUnit = 'kPa',Lower = 0, Symbol = "\Delta P_{fric}");
    82         dPelv(N+1)      as pressure     (Brief = "Elevation Pressure Drop", DisplayUnit = 'kPa',Lower = 0 , Symbol = "\Delta P_{elev}");
    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}");
    86         f(N+1)                  as fricfactor   (Brief = "Friction Factor");
     77       
     78        Liquid          as props                (Brief = "Liquid Properties", Symbol = " ");
     79        Vapour     as props                     (Brief = "Vapour Properties", Symbol = " ");
     80       
     81        Tincr(N+1)      as temperature (Brief = "Temperature Profile", Symbol = "T_{incr}");
     82       
     83        Pdrop                           as pressure             (Brief = "Total Pressure Drop", DisplayUnit = 'kPa',Lower = 0, Symbol = "\Delta P_{drop}");
     84        dPfricTotal             as pressure             (Brief = "Total Friction Pressure Drop", DisplayUnit = 'kPa',Lower = 0, Symbol = "\Delta Ptotal_{fric}");
     85        dPelvTotal              as pressure             (Brief = "Total Elevation Pressure Drop", DisplayUnit = 'kPa',Lower = 0 , Symbol = "\Delta Ptotal_{elev}");     
     86        dPaccTotal              as pressure             (Brief = "Total Acceleration Pressure Drop", DisplayUnit = 'kPa',Lower = 0 , Symbol = "\Delta Ptotal_{acc}");   
     87        dPfric(N+1)             as pressure             (Brief = "Friction Pressure Drop", DisplayUnit = 'kPa',Lower = 0, Symbol = "\Delta P_{fric}");
     88        dPelv(N+1)              as pressure             (Brief = "Elevation Pressure Drop", DisplayUnit = 'kPa',Lower = 0 , Symbol = "\Delta P_{elev}");
     89        dPacc(N+1)              as pressure             (Brief = "Acceleration Pressure Drop", DisplayUnit = 'kPa',Lower = 0 , Symbol = "\Delta P_{acc}");
     90        Pincr(N+1)              as pressure             (Brief = "Pressure Profile", DisplayUnit = 'kPa' , Symbol = "P_{incr}");
     91       
     92        Lincr(N+1)              as length                       (Brief = "Length Points", Symbol = "L_{incr}");
     93        f(N+1)                          as fricfactor           (Brief = "Friction Factor");
     94        Mw                                      as molweight    (Brief = "Average Mol Weight");
     95        zmass(NComp)    as fraction                     (Brief = "Mass Fraction");
     96        Re(N+1)                         as Real                         (Brief = "Reynolds Number Profile");
     97        Fw                                      as flow_mass    (Brief = "Total Mass Flow Profile" , Lower = 0,DisplayUnit = 'kg/h', Symbol = "F_w");
     98       
     99        rho(N+1)                        as dens_mass    (Brief = "Mass Density Profile" , Symbol = "\rho");
     100        mu(N+1)                         as viscosity            (Brief = "Viscosity Profile" , Symbol = "\mu");
     101        Vel(N+1)                        as velocity             (Brief = "Velocity Profile");
     102        vm(N+1)                 as vol_mol              (Brief = "Mixture Molar Volume Profile");
     103        h(N+1)                          as enth_mol             (Brief = "Molar Enthalpy profile");
     104       
     105#       dPdLfric(N+1)           as positive             (Brief = "Friction Gradient",Lower = 0, Unit = 'Pa/m', DisplayUnit = 'kPa/m');
     106#       dPelvdL(N+1)            as positive             (Brief = "Elevation Gradient", Lower = 0 , Unit = 'Pa/m', DisplayUnit = 'kPa/m');
     107#       dPaccdL(N+1)            as positive             (Brief = "Acceleration Gradient",Lower = 0 , Unit = 'Pa/m', DisplayUnit = 'kPa/m');
     108       
     109#Beggs&Brill Method
     110#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     111        NFR(N+1)                                        as positive             (Brief="Froude Number", Lower=1E-10,Symbol = "N_{FR}");
     112        BBLine1(N+1)                    as positive             (Brief="Beggs and Brill Correlation Parameter 1");
     113        BBLine2(N+1)                    as positive             (Brief="Beggs and Brill Correlation Parameter 2");
     114        BBLine3(N+1)                    as positive             (Brief="Beggs and Brill Correlation Parameter 3");
     115        BBLine4(N+1)                    as positive             (Brief="Beggs and Brill Correlation Parameter 4");
     116        Map(N+1)                                        as positive             (Brief="Beggs and Brill Flag ");#apenas verificação do algoritmo
     117        Atrans(N+1)                             as positive             (Brief="Beggs and Brill Correction When Flow is in Transition Pattern");
     118        Y(N+1)                                          as positive             (Brief="Beggs and Brill Correction for Friction Factorin Two Phase Flow");
    87119       
    88120EQUATIONS
    89121
    90 "Mass Flow"
    91         Properties.Fw = sum(M*Inlet.z)*Inlet.F;
     122"Total Average Molecular Weight"
     123        Mw = sum(M*Inlet.z);
     124
     125"Froude Number"
     126        NFR =Vel*Vel/(g*Dpipe);
     127
     128"No Slip Liquid Volume Fraction"
     129        Liquid.Vsfrac = Liquid.Vel/Vel;
     130
     131"No Slip Vapour Volume Fraction"
     132        Vapour.Vsfrac = Vapour.Vel/Vel;
     133
     134"Vapour Holdup"
     135        Vapour.Holdup=1-Liquid.Holdup;
     136       
     137for i in [1:N+1]
     138
     139if Vapour.frac(i) equal 0 then # Only Liquid phase
     140
     141"Reynolds Number"
     142        Re(i)*mu(i)     = rho(i)*Vel(i)*Dpipe;
     143
     144"Liquid Superficial Velocity"
     145        Liquid.Vel(i) = Liquid.Fw(i)/Apipe/Liquid.rho(i);
     146
     147"Vapour Superficial Velocity"
     148        Vapour.Vel(i) = 0*'m/s';
     149
     150"Liquid Mass Flow"
     151        Liquid.Fw(i) = Fw;
     152
     153"Vapour Mass Flow"
     154        Vapour.Fw(i) =  0*Fw;
     155       
     156"Liquid Holdup"
     157        Liquid.Holdup(i) = 1;
     158
     159"Beggs and Brill Line 1"
     160        BBLine1(i) = 1;
     161
     162"Beggs and Brill Line 2"
     163        BBLine2(i) = 1;
     164
     165"Beggs and Brill Line 3"
     166        BBLine3(i) = 1;
     167
     168"Beggs and Brill Line 4"
     169        BBLine4(i) = 1;
     170       
     171Map(i)=10;     
     172Atrans(i) = 1;
     173Y(i) = 1;
     174
     175else if Liquid.frac(i) equal 0 then # Only Vapour phase
     176
     177"Reynolds Number"
     178        Re(i)*mu(i)     = rho(i)*Vel(i)*Dpipe;
     179
     180"Vapour Superficial Velocity"
     181        Vapour.Vel(i) = Vapour.Fw(i)/Apipe/Vapour.rho(i);
     182
     183"Liquid Superficial Velocity"
     184        Liquid.Vel(i) = 0*'m/s';
     185
     186"Liquid Mass Flow"
     187        Liquid.Fw(i) = 0*Fw;
     188
     189"Vapour Mass Flow"
     190        Vapour.Fw(i) = Fw;
     191
     192"Liquid Holdup"
     193        Liquid.Holdup(i) = 0;
     194
     195"Beggs and Brill Line 1"
     196        BBLine1(i) = 1;
     197
     198"Beggs and Brill Line 2"
     199        BBLine2(i) = 1;
     200
     201"Beggs and Brill Line 3"
     202        BBLine3(i) = 1;
     203
     204"Beggs and Brill Line 4"
     205        BBLine4(i) = 1;
     206
     207Map(i)=11;     
     208Atrans(i) = 1;
     209Y(i) = 1;
     210
     211else # Two phase
     212
     213"Reynolds Number"
     214        Re(i)*(Liquid.mu(i)*Liquid.Vsfrac(i)+Vapour.mu(i)*Vapour.Vsfrac(i))= (Liquid.rho(i)*Liquid.Vsfrac(i)+Vapour.rho(i)*Vapour.Vsfrac(i))*Vel(i)*Dpipe; #no sip Reynolds Number
     215
     216"Vapour Superficial Velocity"
     217        Vapour.Vel(i) = Vapour.Fw(i)/Apipe/Vapour.rho(i);
     218
     219"Liquid Superficial Velocity"
     220        Liquid.Vel(i) = Liquid.Fw(i)/Apipe/Liquid.rho(i);
     221
     222"Liquid Mass Flow"
     223        Liquid.Fw(i)= Inlet.F*Liquid.frac(i)*sumt(M(1:NComp)*Liquid.Mfrac(i,1:NComp));
     224
     225"Vapour Mass Flow"
     226        Vapour.Fw(i) =  Inlet.F*Vapour.frac(i)*sumt(M(1:NComp)*Vapour.Mfrac(i,1:NComp));
     227
     228switch Correlation
     229       
     230        case "Beggs-Brill":
     231       
     232"Beggs and Brill Line 1"       
     233        BBLine1(i) = 2.499687083 + 0.302*log(Liquid.Vsfrac(i));#aplicado o log somente no lado direito da equação
     234
     235"Beggs and Brill Line 2"       
     236        BBLine2(i) = -3.033764376 - 2.4684*log(Liquid.Vsfrac(i));#aplicado o log somente no lado direito da equação
     237
     238"Beggs and Brill Line 3"       
     239        BBLine3(i) = -1 - 1.4516*log(Liquid.Vsfrac(i));#aplicado o log somente no lado direito da equação
     240
     241"Beggs and Brill Line 4"       
     242        BBLine4(i) = -0.301029996 - 6.738*log(Liquid.Vsfrac(i));#aplicado o log somente no lado direito da equação
     243
     244Y(i)*Liquid.Holdup(i)*Liquid.Holdup(i) = Liquid.Vsfrac(i);
     245
     246# Find the Flow Pattern for Beggs-Brill Method
     247#++++++++++++++++++++++++++++++++++++++++++
     248        if NFR(i) < 10^BBLine1(i) and Liquid.Vsfrac(i) < 0.01
     249        then
     250        "Liquid Holdup in Segregated Flow"
     251        Liquid.Holdup(i)*NFR(i)^0.0868 = 0.98*Liquid.Vsfrac(i)^(0.4846);
     252        Atrans(i) = 1;
     253        Map(i)=1;
     254
     255        else
     256        if  NFR(i) < 10^BBLine2(i) and Liquid.Vsfrac(i) >= 0.01
     257        then
     258        "Liquid Holdup in Segregated Flow"
     259        Liquid.Holdup(i)*NFR(i)^0.0868 = 0.98*Liquid.Vsfrac(i)^(0.4846);
     260        Atrans(i) = 1;
     261        Map(i)=1;
     262
     263        else
     264        if NFR(i) > 10^BBLine2(i) and NFR(i) <= 10^BBLine3(i) and Liquid.Vsfrac(i) >= 0.01
     265        then
     266        "Liquid Holdup in Transition Flow"#usando segregated, deve ser por interpolação!!!
     267        Liquid.Holdup(i)*NFR(i)^0.0868 = 0.98*Liquid.Vsfrac(i)^(0.4846);
     268       
     269        Atrans(i) = abs(10^BBLine3(i) - NFR(i))/(abs(10^BBLine3(i) - 10^BBLine2(i))+eps);
     270        Map(i)=2;
     271
     272        else
     273        if NFR(i) > 10^BBLine3(i) and NFR(i) <= 10^BBLine1(i) and Liquid.Vsfrac(i) >= 0.01 and Liquid.Vsfrac(i) < 0.4
     274        then
     275        "Liquid Holdup in Intermittent Flow"
     276        Liquid.Holdup(i)*NFR(i)^0.0173 = 0.845*Liquid.Vsfrac(i)^(0.5351);
     277        Atrans(i) = 1;
     278        Map(i)=3;
     279
     280        else
     281        if NFR(i) > 10^BBLine3(i) and NFR(i) <= 10^BBLine4(i) and Liquid.Vsfrac(i) >=  0.4
     282        then
     283        "Liquid Holdup in Intermittent Flow"
     284        Liquid.Holdup(i)*NFR(i)^0.0173 = 0.845*Liquid.Vsfrac(i)^(0.5351);
     285        Atrans(i) = 1;
     286        Map(i)=3;
     287
     288        else
     289        if NFR(i) >= 10^BBLine1(i) and Liquid.Vsfrac(i) <= 0.4
     290        then
     291        "Liquid Holdup in Distributed Flow"
     292        Liquid.Holdup(i)*NFR(i)^0.0609 = 1.065*Liquid.Vsfrac(i)^(0.5824);
     293        Atrans(i) = 1;
     294        Map(i)=4;
     295
     296        else
     297        if NFR(i) > 10^BBLine4(i) and Liquid.Vsfrac(i) >= 0.4
     298        then
     299        "Liquid Holdup in Distributed Flow"
     300        Liquid.Holdup(i)*NFR(i)^0.0609 = 1.065*Liquid.Vsfrac(i)^(0.5824);
     301        Atrans(i) = 1;
     302        Map(i)=4;
     303
     304        else
     305        "Outside the Range of Beggs-Brill Method"
     306        Liquid.Holdup(i) = Liquid.Vsfrac(i);
     307        Atrans(i) = 1;
     308        Map(i)=5;
     309
     310end
     311
     312end
     313
     314end
     315
     316end
     317
     318end
     319
     320end
     321
     322end
     323#++++++++++++++++++++++++++++++++++++++++++
     324        case "Lockhart-Martinelli":
     325       
     326        BBLine1(i) = 1;
     327
     328        BBLine2 (i)= 1;
     329
     330        BBLine3(i) = 1;
     331
     332        BBLine4 (i)= 1;
     333
     334"Liquid Holdup"
     335        Liquid.Holdup(i)*NFR(i)^0.0868 = 0.98*Liquid.Vsfrac(i)^(0.4846);#usando beggs-brill
     336        Atrans(i) = 1;
     337        Map(i)=20;
     338        Y(i)=1;
     339
     340end
     341
     342
     343end
     344
     345end
     346
     347end
     348
     349"Mass Fraction"
     350        zmass = M*Outlet.z / Mw;       
     351
     352"Total Mass Flow"
     353        Fw = sum(M*Inlet.z)*Inlet.F;
    92354
    93355"Inlet Boudary for Temperature Profile"
     
    96358"Outlet Boundary for Temperature Profile"
    97359        Tincr(N+1) = Outlet.T;
    98        
     360
    99361"Inlet Boudary for Pressure Profile"
    100362        Pincr(1) =  Inlet.P;
     
    120382"Outlet Composition"
    121383        Outlet.z = Inlet.z;
    122        
     384
    123385"Molar Balance"
    124386        Outlet.F = Inlet.F;
    125387
    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);
     388"Mixture Velocity"
     389        Vel= Vapour.Vel+Liquid.Vel;
    131390
    132391"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);
     392        #dPfric(1:N+1)  = (0.5*f(1:N+1)*Lincr(1:N+1)*rho(1:N+1)*Vel(1:N+1)^2/Dpipe);
     393        dPfric  = (0.5*f*Lincr*(Liquid.rho*Liquid.Vsfrac+Vapour.rho*Vapour.Vsfrac)*Vel*Vel/Dpipe);
     394       
     395        #DPDLfric       = (0.5*f*(Liquid.rho*Liquid.Vsfrac+Vapour.rho*Vapour.Vsfrac)*Vel^2/Dpipe);
    134396
    135397"Incremental Acceleration Pressure Drop at Pipe Entry"
    136398        dPacc(1) = 0*'Pa';
    137399
    138 if 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
    145400"Incremental Elevation Pressure Drop"
    146         dPelv(1:N+1)    = Properties.rho(1:N+1)*g*Lincr(1:N+1)*(Hrise/abs(Lpipe));
    147        
    148 end
    149 
     401        dPelv= rho*g*Lincr*sin(Hrise);
     402       
    150403for i in [1:N+1]
    151404
    152405"Flash Calculation"
    153         [Properties.Vfrac(i), Properties.x(i,:), Properties.y(i,:)] = PP.FlashPH(Pincr(i), Properties.h(i), Inlet.z);
     406        [Vapour.frac(i), Liquid.Mfrac(i,:), Vapour.Mfrac(i,:)] = PP.FlashPH(Pincr(i), h(i), Inlet.z);
    154407
    155408"Liquid Fraction"
    156         Properties.Lfrac(i) = 1-Properties.Vfrac(i);
     409        Liquid.frac(i) = 1-Vapour.frac(i);
    157410
    158411"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,:));
     412        h(i) = Liquid.frac(i)*Liquid.h(i) + Vapour.frac(i)*Vapour.h(i);
     413       
     414        Liquid.h(i) = PP.LiquidEnthalpy(Tincr(i),Pincr(i),Liquid.Mfrac(i,:));
     415       
     416        Vapour.h(i) = PP.VapourEnthalpy(Tincr(i),Pincr(i),Vapour.Mfrac(i,:));
    160417
    161418"Density"
    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,:));
    163        
    164 "Viscosiyty"
    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,:));
     419        rho(i) = Liquid.frac(i)*Liquid.rho(i)+Vapour.frac(i)*Vapour.rho(i);
     420       
     421        Liquid.rho(i) = PP.LiquidDensity(Tincr(i),Pincr(i),Liquid.Mfrac(i,:));
     422       
     423        Vapour.rho(i) = PP.VapourDensity(Tincr(i),Pincr(i),Vapour.Mfrac(i,:));
     424       
     425"Viscosiyty"
     426       
     427        mu(i) = Liquid.frac(i)*Liquid.mu(i)+Vapour.frac(i)*Vapour.mu(i);
     428       
     429        Liquid.mu(i) = PP.LiquidViscosity(Tincr(i),Pincr(i),Liquid.Mfrac(i,:));
     430       
     431        Vapour.mu(i) = PP.VapourViscosity(Tincr(i),Pincr(i),Vapour.Mfrac(i,:));
    166432       
    167433"Molar Volume"
    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,:));
     434        vm(i) = Liquid.frac(i)*Liquid.vm(i)+Vapour.frac(i)*Vapour.vm(i);
     435       
     436        Liquid.vm(i) = PP.LiquidVolume(Tincr(i),Pincr(i),Liquid.Mfrac(i,:));
     437       
     438        Vapour.vm(i) = PP.VapourVolume(Tincr(i),Pincr(i),Vapour.Mfrac(i,:));
    169439
    170440end
     
    176446
    177447"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);
     448        #dPacc(i+1) = (Fw/Apipe)^2*(1/(Liquid.rho(i+1)*Liquid.Holdup(i+1)+Vapour.rho(i+1)*Vapour.Holdup(i+1))-1/
     449        #(Liquid.rho(i)*Liquid.Holdup(i)+Vapour.rho(i)*Vapour.Holdup(i)));
     450        dPacc(i+1) = (Fw/Apipe)^2*(1/rho(i+1)-1/rho(i));
     451        #dPacc(i+1) = dPacc(1);
    180452       
    181453"Incremental Length"
     
    189461
    190462        case "Constant Temperature":
    191                 "Outlet Temperature"
     463        "Outlet Temperature"
    192464                Tincr(i+1) = Inlet.T;
    193465       
    194466        case "Linear Temperature profile":
    195                 "Outlet Temperature"
     467        "Outlet Temperature"
    196468                Toutlet - Tincr(1) = Lpipe*((Tincr(i+1)-Tincr(i))/(Lincr(i+1)-Lincr(i)));
    197469
     
    201473
    202474for i in [1:N+1]
    203        
     475
    204476switch FlowRegime
    205        
     477
    206478        case "laminar":
    207        
     479
    208480"Friction Factor for Pressure Drop - laminar Flow"
    209         f(i)*Properties.Re(i) = 16;
    210        
    211         when Properties.Re(i) > 2300 switchto "turbulent";
     481        f(i)*Re(i) = 16;
     482
     483        when Re(i) > 2300 switchto "turbulent";
    212484
    213485        case "turbulent":
    214486
    215487"Friction Factor  for Pressure Drop - Turbulent Flow"
    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        
    220 end
    221 
    222 end
    223 
    224 end
     488        1/sqrt(f(i))= -2*log(Roughness/Dpipe/3.7+2.51/Re(i)/sqrt(f(i)));
     489
     490        when Re(i) <= 2300 switchto "laminar";
     491
     492end
     493
     494end
     495
     496end
     497
     498FlowSheet Pipe_simples
     499
     500PARAMETERS
     501        PP      as Plugin(Brief="Physical Properties",
     502                Type="PP",
     503                Components = [ "isobutane", "benzene","n-hexane"],
     504                LiquidModel     = "PR",
     505                VapourModel     = "PR"
     506        );
     507
     508        NComp   as Integer;
     509
     510DEVICES
     511        Tube as pipe;
     512        Feed as simple_source;
     513        Storage as simple_sink;
     514
     515SET
     516        NComp  = PP.NumberOfComponents;
     517        Tube.N  = 2;
     518       
     519        #Tube.Thermal ="Constant Temperature";
     520        Tube.Thermal ="Linear Temperature profile";
     521        Tube.Toutlet  = (80+273.15) * 'K';
     522       
     523       
     524CONNECTIONS
     525        Feed.Outlet to Tube.Inlet;
     526        Tube.Outlet to Storage.Inlet;
     527
     528SPECIFY
     529        Feed.Outlet.F = 10 * 'kmol/h';
     530        Feed.Outlet.P = 3* 'atm';
     531        Feed.Outlet.T = (30+273.15) * 'K';
     532        Feed.Outlet.z = [0.5,0.2,0.3];
     533
     534SET
     535        Tube.Lpipe = 5*'m';
     536        Tube.Hrise =    0*'deg'; # working only in the horizontal , must be updated
     537        Tube.Correlation = "Beggs-Brill";
     538       
     539
     540        Tube.Dpipe                      = 3*'in';
     541        Tube.Roughness  = 4.572E-05*'m';
     542       
     543        OPTIONS
     544        Dynamic = false;
     545        #GuessFile="Pipe_simples";
     546
     547end
Note: See TracChangeset for help on using the changeset viewer.