Transcript A Stability/Bifurcation Framework For Process Design C. Theodoropoulos1, N. Bozinis2, C. Siettos1, C.C.

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.
parameter
• RPM: method of choice to build around existing time-stepping codes.
–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)
Reconstruct solution:
u = p+q = PN(p,q)+QF
Initial state un
F.P.I.
Newton
iterations
Timestepping
Legacy Code
F(un)
Subspace
Q =I-P
Subspace P of few NO
slow eigenmodes
• Treats timstepping routine, as a
“black-box”
– Timestepper evaluates un+1= F(un)
•
Recursively identifies subspace of
slow eigenmodes, P
•
Substitutes pure Picard iteration with
–Newton method in P
–Picard iteration in Q = I-P
Picard
iteration
Convergence?
YES
Final state uf
•
Reconstructs solution u from sum of
the projectors P and Q onto subspace
P and its orthogonal complement Q,
respectively:
–u = PN(p,q) + QF
Princeton University
gPROMS:A General Purpose Package
Nonlinear
algebraic
equation
solvers
Differential
algebraic
equation
solvers
gPROMS
Steadystate &
Dynamic
Simulation
Parameter
Estimation
Nonlinear
programming
solvers
Steady-state
& Dynamic
Optimisation
Dynamic
optimisation
solvers
gPROMS
Model
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
•well-posedness
 structural analysis algorithms
•DAE index analysis
•consistency of DAE IC’s
•automatic block triangularisation
Advanced features:
 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
2
ut   u  u  u 3  v
– e.g. Beloushov-Zhabotinski
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 Turning Point
<u>
Around Hopf
e
Princeton University
Eigenspectrum Around Hopf
Princeton University
Eigenvectors
e = 0.02
Princeton University
Arc-length continuation with gPROMS
System:
Det[
0  f ( y, p)
dy
 f ( y, p)
dt
y
f ( y; p*)
]0
y
p
Pseudo – arc length condition
(1)
( y1  y0 )T
( p  p0 )
( y  y1 )  1
( p  p1 )  S  0
S
S
(2)
Solve (1) & (2)
continuation
(II)
through
FORTRAN
F.P.I
continuation
(I)
within
gPROMS
Princeton University
System Jacobian
ODEs :
dx
 f ( x, p )
dt
dx
 f ( x, y , p )
DAEs : dt
y  y * ( x)
0  g ( x, y , p )
F.P.I
Continuation
within
gPROMS
F.P.I
R.P.M.
through
FORTRAN
Getting system
Jacobian
through an FPI
Obtain “correct”
Jacobian of leading
eigenspectrum
Cannot get “correct”
Jacobian from
augmented system
1
f ( x , p )
x
f f  g  g
  
x y  y  x
Jacobian of the ODE
Stability matrix
Princeton University
Tubular Reactor: A DAE system
Dimensionless equations:
2
x1
x1
x2
1  x1
 Pe1


Da
(
1

x
)
exp[
]
1
2
t
z
z
1  x2 / 
(1)
2
x2
x2
x2
1  x2
 Pe2

 x2  BDa(1  x1 ) exp[
]  x2 w
2
t
z
z
1  x2 / 
(2)
Boundary Conditions:
x1 ( z  0, t )
 Pe1 x1  0
z
x1 ( z  1, t )
0
z
x2 ( z  0, t )
 Pe 2 x2  0
z
x2 ( z  1, t )
0
z
(3)
(4)
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
•101-node FD discretization
1
•2 unknowns (x1,x2) per node
0.8
Hopf pt.
•State variables:
99 (x 2) unknowns at inner nodes
•Perform RPM-gPROMs at 99-space
to obtain correct Jacobian
x1
0.6
0.4
0.2
0
0.1
0.11
0.12
0.13
0.14
Da
Princeton University
Eigenspectrum
1
6.00E-02
5.00E-02
0.5

4.00E-02
0
-1
-0.5
0
0.5
1
1.5
3.00E-02
-0.5
2.00E-02
-1
1.00E-02

0.00E+00
Da=0.121738
0
1
20
40
60
80
100
120
3.50E-02
3.00E-02
0.5
2.50E-02
2.00E-02
Im
0
-1
-0.5
0
0.5
-0.5
1
1.5
1.50E-02
Re
1.00E-02
5.00E-03
-1
0.00E+00
0
20
40
60
80
100
120
Da=0.110021
Princeton University
Stability Analysis without the Equations
SYSTEM AROUND STEADY STATE
y k 1  (yk )
1
Leading
Spectrum
0.5
0
-1
-0.5
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
Aq 
(y k  e q)  (y k )
e
(2) Compute and store ht , j   Aq j , q t  , t  1, 2,... j
j
(3) r j  Aq j   ht , j q t
t 1
Large-scale eigenvalue calculations
(Arnoldi using system Jacobian):
R.B. Lechouq & A.G. Salinger,
Int. J. Numer. Meth.(2001)
(4) h j 1, j   r j , r j 
q  r j / h j 1, j
(5) j 1
1/ 2
End For
Princeton University
Rapid Pressure Swing Adsorption
1-Bed 2-Step Periodic Adsorption Process
t=0 to T/2
Ci(z=0)=PfYf/(RTf)
P(z=0)=Pf
z=L
z=0
•Isothermal operation
•Modeling Equations (Nilchan & Pantelides)
Mass balance in ads. bed
Ci
qi
 (vCi )
 2Ci
et
 b

 Di
t
t
z
z 2
n
P
  Ci
RT i 1
P
180v (1  e b ) 2

z
d p2
e b3
qi
 ki (mi pi  qi )
t
Step 1 :
Pressurisation
t= T/2 to T
Ci
( z  0)  0
z
P(z=0)=Pw
Darcy’s law
Rate of ads
.
Step 2:
Depressurisation
Princeton University
Rapid Pressure Swing Adsorption
1-Bed 2-Step Periodic Adsorption Process
q , c (t=T)
Production
Zeolite
Bed
of oxygen enriched air
5A adsorbent (300m)
q ,c (t=0)
1m long, 5cm diameter
Short
cycle
q , c (t=T/2)
Must obtain:
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
0.5
50
c1(z=0.5)
45
(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
0
1
0.2
0.4
z
0.6
0.8
x
1
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
z
0.8
x
1
0
0.2
0.4
0.6
0.8
x
1
z
Princeton University
Leading Eigenvectors, =0.99484
0.0012
0.16
c1
0.12
q1
q1
c1
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
-0.15
q2
q2
c2
-6.00E-04
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
(II)
through
FORTRAN
F.P.I
continuation
(I)
within
gPROMS
F.P.I
F.P.I
tracing
limit cycles
within gPROMS
R.P.M
through
FORTRAN
Getting system
Jacobian
through an FPI
Princeton University
Tracing limit cycles
Tracing Limit Cycles
SYSTEM:
dy
 f ( y, p)
dt
Periodic Solutions:
y(t+T)=y(t)
Period T not known beforehand
d y
 dt   f ( y, p)
 dT   

0




 dt 
y(0)  y(T )  0
G( y(0), p)  0
G( y(0), p)  yi (0)  a  0
dy (0)
G ( y (0), p )  i
0
dt
Princeton University