VIEWS: 22 PAGES: 24 POSTED ON: 12/27/2010
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 180v (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 (300m) 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