Docstoc

Princeton activities

Document Sample
Princeton activities Powered By Docstoc
					A Stability/Bifurcation Framework
       For Process Design

C. Theodoropoulos1, N. Bozinis2, C. Siettos1,
   C.C. Pantelides2 and I.G. Kevrekidis1

         1Department
                   of Chemical Engineering,
     Princeton University, Princeton, NJ 08544

     2 Centre for Process System Engineering,
     Imperial College, London, SW7 2BY, UK

                                            Princeton University
                                   Motivation
• A large number of existing scientific, large-scale legacy codes
   –Based on transient (timestepping) schemes.
• Enable legacy codes perform tasks such as bifurcation/stability analysis
   –Efficiently locate multiple steady states and assess the stability of solution branches.
   –Identify the parametric window of operating conditions
    for optimal performance
   –Locate periodic solutions
     •Autonomous, forced (PSA,RFR)
   –Appropriate controller design.
• RPM: method of choice to build around existing time-stepping codes.               parameter

   –Identifies the low-dimensional unstable subspace of a few “slow” eigenvalues
   –Stabilizes (and speeds-up) convergence of time-steppers even onto unstable steady-states.
   –Efficient bifurcation analysis by computing only the few eigenvalues of the small
   subspace.
     •Even when Jacobians are not explicitly available (!)

                                                                   Princeton University
           Recursive Projection Method (RPM)
                                                       • Treats timstepping routine, as a
 Reconstruct solution:                                   “black-box”
                                  Initial state un
u = p+q = PN(p,q)+QF                                       – Timestepper evaluates un+1= F(un)
                                           F.P.I.      •   Recursively identifies subspace of
                 Newton           Timestepping
                                                           slow eigenmodes, P
                 iterations       Legacy Code          •   Substitutes pure Picard iteration with
                                   F(un)
                                            Picard
                                                            –Newton method in P
                                           iteration        –Picard iteration in Q = I-P
Subspace   Subspace P of few NO
Q =I-P      slow eigenmodes
                                   Convergence?        •   Reconstructs solution u from sum of
                                                           the projectors P and Q onto subspace
                                            YES            P and its orthogonal complement Q,
                                                           respectively:
                                   Final state uf           –u = PN(p,q) + QF


                                                                        Princeton University
               gPROMS:A General Purpose Package
Nonlinear                                                                Nonlinear
algebraic                                                              programming
equation                                               gPROMS             solvers
 solvers            Steady-                           Steady-state
                     state &                           & Dynamic
Differential
                    Dynamic                           Optimisation
algebraic          Simulation
 equation                                                             Dynamic
  solvers
                                  gPROMS                             optimisation
                    Parameter       Model                              solvers
                    Estimation
                                                      Data
                                                  Reconciliation



                                 Maximum likelihood
                                 estimation solvers


                                                            Princeton University
    Mathematical solution methods in gPROMS
•   Combined symbolic, structural & numerical techniques
      symbolic differentiation for partial derivatives
      automatic identification of problem sparsity
      structural analysis algorithms                   •well-posedness
                                                            •DAE index analysis
•   Advanced features:                                      •consistency of DAE IC’s
                                                            •automatic block triangularisation
      exploitation of sparsity at all levels
      support for mixed analytical/numerical partial derivatives
      handling of symmetric/asymmetric discontinuities at all levels
•   Component-based architecture for numerical solvers
      open interface for external solver components
      hierarchical solver architectures
          • mix-and-match
          • external solvers can be introduced at any level of the hierarchy

                                                                      Princeton University
  FitzHugh-Nagumo: An PDE-based Model

• Reaction-diffusion model in one dimension
• Employed to study issues of pattern formation
  in reacting systems
   – e.g. Beloushov-Zhabotinski
                                   ut   u  u  u 3  v
                                       2


     reaction                         vt  δ 2v  ε(u  a1v  a0 )
   – u “activator”, v “inhibitor”
   – Parameters: δ  4.0, a0  0.03, a1  2.0
   – no-flux boundary conditions
   – e, time-scale ratio, continuation parameter
• Variation of e produces turning points
  and Hopf bifurcations
                                                    Princeton University
           Bifurcation Diagrams
      Around Hopf        Around Turning Point
<u>




            e




                                  Princeton University
Eigenspectrum Around Hopf




                     Princeton University
Eigenvectors
    e = 0.02




               Princeton University
      Arc-length continuation with gPROMS
                        dy
         System:            f ( y, p )
                        dt                                                    y
                 f ( y; p*)
         Det [               ] 0
                      y                                                                   p
                                 Pseudo – arc length condition
0  f ( y, p)    (1)      ( y1  y0 )T             ( p  p0 )
                                       ( y  y1 )  1         ( p  p1 )  S  0   (2)
                               S                     S

                              Solve (1) & (2)
                                                    continuation
 continuation                 F.P.I                      (I)
                                                      within
      (II)
                                                     gPROMS
   through
  FORTRAN

                                                                           Princeton University
                        System Jacobian
                                       dx
             dx                            f ( x, y , p )
ODEs :           f ( x, p )    DAEs : dt                  y  y * ( x)
             dt                        0  g ( x, y , p )

                               F.P.I        R.P.M.              Obtain “correct”
                                           through              Jacobian of leading
                                          FORTRAN               eigenspectrum
             Continuation
               within          F.P.I     Getting system          Cannot get “correct”
              gPROMS                        Jacobian             Jacobian from
                                         through an FPI          augmented system

                                                                  1
         f ( x, p )                            f f  g  g
                                                    
            x                                  x y  y  x
                                                       
  Jacobian of the ODE                               Stability matrix

                                                          Princeton University
      Tubular Reactor: A DAE system
Dimensionless equations:
  x1       1  x1
                2
                     x                       x2
       Pe1          1  Da(1  x1 ) exp[           ]                        (1)
  t           z 2
                     z                   1  x2 / 
  x2       1  x2
                2
                      x2                              x2
       Pe2               x2  BDa(1  x1 ) exp[           ]  x2 w       (2)
   t          z 2
                      z                           1  x2 / 
Boundary Conditions:
      x1 ( z  0, t )                   x2 ( z  0, t )
                        Pe1 x1  0                        Pe2 x2  0        (3)
            z                                 z
       x1 ( z  1, t )                  x2 ( z  1, t )
                        0                                0                  (4)
              z                               z
 Eqns (1)-(4): system of DAEs. Can also substitute to obtain system of ODEs.

                                                            Princeton University
           Bifurcation/Stability with RPM-gPROMS

                                                 •Model solved as DAE system
                                                     •2 algebraic equations @ each boundary
      1                                          •101-node FD discretization

     0.8                                         •2 unknowns (x1,x2) per node
                 Hopf pt.
     0.6                                         •State variables:
x1




     0.4
                                                  99 (x 2) unknowns at inner nodes
                                                 •Perform RPM-gPROMs at 99-space
     0.2                                          to obtain correct Jacobian
      0
           0.1     0.11     0.12   0.13   0.14
                            Da




                                                                     Princeton University
                                 Eigenspectrum
              1                         6.00E-02


            0.5                         5.00E-02
                                                                            
                                        4.00E-02
              0
                                        3.00E-02
-1   -0.5          0   0.5   1    1.5                                                  
            -0.5                        2.00E-02


             -1                         1.00E-02

                                        0.00E+00

     Da=0.121738                                    0   20   40        60        80         100    120



              1                          3.50E-02

                                         3.00E-02
            0.5                          2.50E-02

                                         2.00E-02                 Im
              0
                                         1.50E-02
-1   -0.5          0   0.5   1    1.5
                                         1.00E-02                       Re
            -0.5
                                         5.00E-03
             -1                          0.00E+00
                                                    0   20   40        60         80         100    120

      Da=0.110021

                                                                                 Princeton University
              Stability Analysis without the Equations
                                                                                                            1


                 SYSTEM AROUND STEADY STATE                         Leading
                        y k 1  (y k )
                                                                                                          0.5




                                                                    Spectrum                  -1   -0.5
                                                                                                            0
                                                                                                                 0   0.5   1




                                                                                                          -0.5




    y(k)                                                                                                   -1




                                                         Matrix-free ARNOLDI
                                                  Choose q1 with q1  1
+
                                                        For j =1 Until Convergence DO
                                                        (1) Compute and store Aq j
    εq
                         (y k  e q)  (y k )         (2) Compute and store ht , j   Aq j , qt  , t  1,2,... j
                  Aq 
                                    e                                           j

                                                        (3) r j  Aq j   ht , j q t
                                                                               t 1


                                                        (4) h j 1, j   r j , r j 
                                                                                      1/ 2
         Large-scale eigenvalue calculations
         (Arnoldi using system Jacobian):
                                                        (5) q j 1  r j / h j 1, j
         R.B. Lechouq & A.G. Salinger,
         Int. J. Numer. Meth.(2001)                     End For


                                                                                      Princeton University
                    Rapid Pressure Swing Adsorption
                1-Bed 2-Step Periodic Adsorption Process
t=0 to T/2                    •Isothermal operation                                       t= T/2 to T
                              •Modeling Equations (Nilchan & Pantelides)
Ci(z=0)=PfYf/(RTf)                                                                       Ci
                              Mass balance in ads. bed                                       ( z  0)  0
               z=L                                                                        z
P(z=0)=Pf                        Ci      qi     (vCi )       2Ci
                              et      b                Di                            P(z=0)=Pw
                                  t      t        z         z 2
                                     n
                               P
                                     Ci
                              RT i 1
                              P    180v (1  e b ) 2
                                                            Darcy’s law
               z=0            z       2
                                      dp     eb3


                              qi
                                   ki (mi pi  qi )     Rate of ads   .
                              t
                                                                               Step 2:
                Step 1 :
                                                                           Depressurisation
             Pressurisation
                                                                            Princeton University
            Rapid Pressure Swing Adsorption
         1-Bed 2-Step Periodic Adsorption Process
                                          q , c (t=T)
Production      of oxygen enriched air
Zeolite   5A adsorbent (300m)           q ,c (t=0)          q , c (t=T/2)

Bed   1m long, 5cm diameter
                                                Must obtain:
Short   cycle                                   q , c (t=T) = q , c (t=0)

   –1.5s pressurisation, 1.5s depressurisation

   – T= 3s

Low   feed pressure (Pf = 3 bar)
Periodic   steady-state operation
   –reached after several thousand cycles

                                                         Princeton University
                     Typical RPSA simulation results
              (Nilchan and Pantelides, Adsorption, 4, 113-147, 1998)




                                                         1




50                                                     0.5


45
              c1(z=0.5)   (mol/m3)
40
                                                         0
                                           -1   -0.5          0   0.5   1
35


30
                                                       -0.5

25
              Time (s)
20
     0   50    100        150        200                -1




                                                                            Princeton University
         PRM-gPROMS Spatial Profiles (t=T)
30
                                               0.3

               c1 mol/m3                                       q1 mol/kg
20                                             0.2




10                                             0.1




                                                0
0
                                                     0   0.2   0.4         0.6   0.8   x   1
     0   0.2   0.4         0.6   0.8   x   1

                     z                                               z
90                                             0.3
               c2 mol/m3                                         q2 mol/kg

60                                             0.2




30                                             0.1




0                                               0
     0   0.2   0.4         0.6   0.8   x   1         0   0.2   0.4         0.6   0.8   x   1

                     z                                               z

                                                                     Princeton University
                  Leading Eigenvectors, =0.99484
0.16                               0.0012

        c1                                     q1
0.12         c1                                q1
                                   0.0008


0.08

                                   0.0004
0.04



   0                                    0



   0                               0.00E+00




-0.05                              -2.00E-04




 -0.1                              -4.00E-04


             c2                                 q2
        c2                                      q2
                                   -6.00E-04
-0.15



                                                     Princeton University
                             Conclusions
•   Can construct a RPM-based computational framework around large-scale
    timestepping legacy codes to enable them converge to unstable steady states
    and efficiently perform bifurcation/stability analysis tasks.
      – gPROMS was employed as a really good simulation tool
      – communication with wrapper routines through F.P.I.
          • Both for PDE and DAE-based systems.
•   Have “brought to light” features of gPROMS for continuation around turning
    points and information on the Jacobian and/or stability matrix at steady states
    of systems.
•   Employed matrix-free Arnoldi algorithms to perform stability analysis of
    steady state solutions without having either the Jacobian or even the equations!
•   Used the RPM-based superstructure to speed-up convergence and perform
    stability analysis of an almost singular periodically-forced system
•   Have enabled gPROMS to trace autonomous limit cycles
•   Newton-Picard computational superstructure for autonomous limit cycles.

                                                                Princeton University
                                  gPROMS
•   General purpose commercial package for modelling, optimization and control
    of process systems.
•   Allows the direct mathematical description of distributed unit operations
•   Operating procedures can be modelled
     – Each comprising of a number of steps
        • In sequence, in parallel, iteratively or conditionally .
•   Complex processes: combination of distributed and lumped unit operations
     – Systems of integral, partial differential, ordinary differential and algebraic
       equations (IPDAEs).
     – gPROMS solves using method of lines family of numerical methods.
•   Reduces IPDAES to systems of DAEs.
     – Time-stepping or pseudo-timestepping.
•   Jacobians NOT explicitly available.
     – Cannot perform systematic bifurcation/stability analysis studies.
                                                                 Princeton University
                       Tracing limit cycles
                 Tracing Limit Cycles

                         continuation
                              (I)             F.P.I        R.P.M
                           within                         through
                          gPROMS                         FORTRAN
continuation   F.P.I
     (II)
  through
 FORTRAN                                      F.P.I
                                                         Getting system
                              tracing                       Jacobian
                            limit cycles                 through an FPI
                          within gPROMS




                                                 Princeton University
                      Tracing limit cycles
                 Tracing Limit Cycles
              dy
 SYSTEM:          f ( y, p )
              dt

Periodic Solutions:   y(t+T)=y(t)

Period T not known beforehand




        d y
         dt   f ( y, p )        y(0)  y (T )  0
         dT             
                0               G( y(0), p)  0
         dt 
                                           G ( y (0), p)  yi (0)  a  0
                                                             dyi (0)
                                           G ( y (0), p)            0
                                                               dt
                                                        Princeton University

				
DOCUMENT INFO
Shared By:
Categories:
Stats:
views:22
posted:12/27/2010
language:English
pages:24