CSTR: Multiplicidade de Soluções e Análise de Estabilidade Argimiro R. Secchi Programa de...

CSTR: Multiplicidade de Soluções e Análise de Estabilidade

Argimiro R. SecchiPrograma de Engenharia Química – COPPE/UFRJ

Rio de Janeiro, RJ

EQE038 – Simulação e Otimização de Processos Químicos

EQ/UFRJmarço de 2014

2

-Multiplicity of steady states

-Linearization

-System stability

-Complex dynamic behaviors (limit cycles, strange attractors)

-Parametric sensitivity and input sensitivity

System AnalysisSystem Analysis

3

Non-isothermal CSTR

Fe , CAf , CBf , Tf

Fwe , Twe

Fws , Tw

Fs , CA , CB , T

V , T

A k B

Multiplicity of Steady StatesMultiplicity of Steady States

4

In a non-isothermal continuous stirred tank reactor, with diameter of 3.2 m

and level control, pure reactant is fed at 300 K and 3.5 m3/h with concentration

of 300 kmol/m3. A first order reaction occur in the reactor, with frequency factor

of 89 s-1 and activation energy of 6 x 104 kJ/kmol, releasing 7000 kJ/kmol of

reaction heat. The reactor has a jacket to control the reactor temperature, with

constant overall heat transfer coefficient of 300 kJ/(h.m2.K). Assume constant

density of 1000 kg/m3 and constant specific heat of 4 kJ/(kg.K) in the reaction

medium. The fully-open output linear valve has a constant of 2.7 m2.5/h.

Process description

Multiplicity of Steady StatesMultiplicity of Steady States

5

•perfect mixture in the reactor and jacket;

•negligible shaft work;

•(-rA) = k CA;

•constant density;

•constant overall heat transfer coefficient;

•constant specific heat;

•incompressible fluids;

•negligible heat loss to surroundings;

(internal energy) (enthalpy);

•negligible variation of potential and kinetic energies;

•constant volume in the jacket;

•thin metallic wall with negligible heat capacity.

Model assumptions

Multiplicity of Steady StatesMultiplicity of Steady States

6

dt

dVFF

dt

Vdsef

)(

se FFdt

dV

)( AAsfAeA

AA rVCFCFdt

dVC

dt

dCV

dt

VCd

Mass balance in the reactorOverall:

Component:

(2)VrCCFdt

dCV AAfAe

A )()(

(1)

eF

V (3)

CSTR modeling

Multiplicity of Steady StatesMultiplicity of Steady States

7

(4)

2 2

ˆˆ ˆ ˆ ˆ ˆ ˆ2 2f s

e f f f f s s r s

v vdV U K F U P V gz F U PV gz q q w

dt

ˆ ˆ ˆH U PV

ˆ ˆ( ) ˆ ˆ ˆe f s r

d VH dH dVV H F H F H q q

dt dt dt

Energy balance in the reactor:

where

ˆˆ ˆ( )e f r

dHV F H H q q

dt

qqTTCpFdt

dTVCp rfe )(

CSTR modeling

Multiplicity of Steady StatesMultiplicity of Steady States

8

(5)where

qr = (-Hr) V (-rA)

k = k0 exp(–E/RT)

(-rA) = k CA

V = A h

Fs = x Cv h

Tw = f(T)

(7)

(6)

(8)

(10)

(12)

Temperature control (14)

At = A + D h (11)

q = U At (T – Tw)

A = D2/4 (9)

x = f(h) Level control (13)

CSTR modeling

Multiplicity of Steady StatesMultiplicity of Steady States

9

variable units of measurement

Fe, Fs m3 s-1

V m3

t, s

CA, CAf kmol m-3

rA kmol m-3 s-1

kg m-3

Cp kJ kg-1 K-1

T, Tf, Tw K

qr, q kJ s-1

U kJ m-2 K-1 s-1

At, A m2

h, D m

Cv m2.5 h-1

x –

Hr, E kJ kmol-1

R kJ kmol-1 K-1

k, k0 s-1

Consistency analysis

Multiplicity of Steady StatesMultiplicity of Steady States

10

variables: Fe, Fs, V, t, CA, CAf, rA, , Cp, T, Tf, Tw, qr, q, U, At, A, h, D, Cv, x, Hr, E, R, k, k0, 27constants: , Cp, U, D, Cv, Hr, E, R, k0 9

specifications: t 1driving forces: Fe, Tf, CAf 3

unknown variables: Fs, V, CA, rA, T, Tw, qr, q, A, At, h, x, k, 14

equations: 14

Degree of Freedom = variables – constants – specifications – driving forces – equations = unknown variables – equations = 27 – 9 – 1 – 3 – 14 = 0

Dynamic Degree of Freedom (index < 2) = differential equations = 3

Needs 3 initial condition: h(0), CA(0), T(0) 3

Consistency analysis

Multiplicity of Steady StatesMultiplicity of Steady States

11

Running EMSO

Open MSO file

Multiplicity of Steady StatesMultiplicity of Steady States

12

Consistency Analysis

Results

13

14

The CSTR example at the steady state satisfy:

0

0

( )1( ) ( )

1

E

RTr Aft

f w EP RT

P

H k e CU AT T T T

V CC k e

e

V

F

01

AfA E

RT

CC

k e

Multiplicity of Steady StatesMultiplicity of Steady States

15

Rewriting the energy balance:

( ) ( )R GQ T Q T0

0

( )( )

1

E

RTr Af

G E

RTP

H k e CQ T

C k e

Q T a T bR ( )

aU A

V Ct

P

1

f t w

P

T U A Tb

V C

T

Q

QG

QR3QR2QR1

1 2

3

4

5

GR dQdQ

dT dT

stable:

GR dQdQ

dT dT

unstable:

Multiplicity of Steady StatesMultiplicity of Steady States

16

( , )dx

F t xdt ( ) 0F x

1( 1) ( ) ( ) ( )( ) ( ) , 0,1,2,k k m kx x J x F x k Newton-Raphson:

( )( ) ( )

( )k

k iij

j

F xJ x

x

and 0 1m k

Path FollowingPath Following

Homotopic Continuation: ( ; ) (1 ) ( ) ( ) 0 , 0 1H x p p F x pG x p

(0) (0)( ) ( ) ( )G x J x x x (0)( ) ( ) ( )G x F x F x

affine homotopy

Newton homotopy

Multiples solutions can be obtained by continuously varying the parameter p

Multiplicity of Steady StatesMultiplicity of Steady States

17

Parametric Continuation: [ ( ); ( )] 0F x s p s where s is some parameterization, e.g., path arc length

( ) ( ) 0 , eF F dx dp

x s p s x px p ds ds

F FDF

x p

Frechet derivative

a point (xo, po) is:

( , )o oF x p

x

- Regular if is non-singular

( , )o oF x p

x

- Turning point if is singular and DF has rank = n

( , )o oF x p

x

- Bifurcation if is singular and DF has rank < n

reparameterization

Path FollowingPath FollowingMultiplicity of Steady StatesMultiplicity of Steady States

18

Solutions: 1) CA = 13,13 kmol/m3 and T = 659,46 K 2) CA = 132,87 kmol/m3 and T = 523,01 K 3) CA = 299,86 kmol/m3 and T = 332,72 K

Example: a) execute flowsheet in file CSTR_noniso.mso with initial condition of 578 K and compare with result changing the initial condition to 579 K; b) find the three steady states using file CSTR_sea.mso by changing the initial guess for T and CA (use the section GUESS).

Multiplicity of Steady StatesMultiplicity of Steady States

19

LinearizationLinearization

Generate linearized model at given operating point.

Implicit DAE:

Considering the specification as input, u(t), (SPECIFY section in EMSO):

And identifying the algebraic variables as y(t):

( ,́ , , , ) 0F x x y u t

( ,́ , ) 0F x x t

ˆ ˆ( ,́ , , ) 0F x x u t

20

Differentiating F:

and extracting:

The partition:

Define the linearized system:

Linearization

´ 0x x y uF dx F dx F dy F du

1

x y x u

dx dxF F F F

dy du

1

x y x u

A BF F F F

C D

'x A x Bu

y C x Du

(index < 2)

23

Example: execute the flowsheet in file CSTR_linearize.mso with the option Linearize = true and evaluate the characteristic values of the Jacobian matrix (matrix A). Repeat the example with the value of Cp 10 times smaller, i.e., 0.4 kJ / (kg K). Compare the ratio between the greater and the smaller characteristic values in module.

Non-isothermal CSTR: linearizationNon-isothermal CSTR: linearization

24

Stability AnalysisStability Analysis

)(tx

( )dx

F xdt 0 0( ) ( )x t y t ( ) ( )x t y t

Liapunov Stability: is stable (or Liapunov stable) if, given > 0, there exists a = () > 0, such that, for any other solution, y(t), of

satisfying , then for t > t0.

25

Stability AnalysisStability Analysis

)(tx

0 0( ) ( )x t y t b Asymptotic Stability: is asymptotic stable if Liapunov stable and there exists a constant b > 0 such that, if then 0)()(lim

tytx

t

( ) ( ) ( )y t x t x t Defining deviation variables:

2[ ( )]( ) ( ( ))

dx F x tF x F x t y O y

dt x

( ( ) )dx dx dy

F x t ydt dt dt

Expanding in Taylor series:

Linearization: [ ( )] ( )dy

J x t y A t ydt

26

( )x tFor an equilibrium point = x*, the stability is characterized by the characteristics values of the Jacobian matrix J(x*) = A:

x* is a hyperbolic point if none characteristics values of J(x*) has zero real part.

x* is a center if the characteristics values are pure imaginary. Fixed point non-hyperbolic.

x* is a saddle point, unstable, if some characteristics values have real part > 0 and the remaining have real part < 0.

x* is stable or attractor or sink point if all characteristics values have real part < 0.

x* is unstable or repulsive or source point if at least one characteristic value have real part > 0.

Stability AnalysisStability Analysis

27

For a second-order linear system:2det( ) ( ) det( ) 0A I tr A A

2( ) ( ) 4det( )

2

tr A tr A A

Stability AnalysisStability Analysis

28

0 0( ) , (0)E

RTeAA A A A Af

FdCC C k e C C C

dt V

00

( ) ( )( ) , (0)

E

RTe r A t w

fp p

F H k e C UA T TdTT T T T

dt V C V C

Considering the CSTR example with constant volume:

0 02

0 02

( , )( ) ( )

E E

RT RTeA

E EART RT

r e r A t

p p p

F Ek e k e C

V RTJ C T

H k e F H k e C UAE

C V RT C V C

Stability AnalysisStability Analysis

29

3 4

3 4

1.6458 10 3.4282 10(13.13, 659.46)

2.7542 10 4.8934 10J

4 4

4 4

1.6260 10 3.1753 10(132.87, 523.01)

1.5852 10 4.4509 10J

5 7

8 4

7.2050 10 6.6220 10(299.86, 332.72)

5.9285 10 1.0944 10J

5

4

7.2051 10

1.0944 10

5

4

6.3659 10

3.4614 10

3

4

1.0205 10

1.3604 10

2) Saddle Point, unstable

1) Stable node

3) Stable Node

Stability AnalysisStability Analysis

300 50 100 150 200 250 300 350

300

350

400

450

500

550

600

650

700

750

800

CA

T

file: CSTR_nla/traj_cstr.m

Stability AnalysisStability Analysis

31

Complex Dynamic BehaviorComplex Dynamic Behavior

Tw

Tw

CA

THopf point

Tw = 200,37 K

unstable solutions

stable solutions

CSTR example:

32

t (h)

t (h)

unstable limit cycle

A limit cycle is stable if all characteristics values of exp(J p) (Floquet multipliers) are inside the unitary cycle, where J is the Jacobian matrix in the cycle, p = 2 / is the oscillation period and = |Hopf|.

file: CSTR_auto/cstr_bif.mso

Complex Dynamic BehaviorComplex Dynamic Behavior

33

Interface EMSO-AUTOInterface EMSO-AUTO

parameters

Equation system

Jacobian matrix

First steady-state solution

34

Interface EMSO-AUTOInterface EMSO-AUTO

21 ( )1 11 x tdx t

x t p x t edt

222 13 14 1 x tdx t

x t p x t edt

2 2

2 2

1

1

1 (1 )( )

14 3 14 (1 )

x x

x x

p e p x eJ x

p e p x e

p = 0: x* = (0, 0) (J) = (-1, -3)

35

Interface EMSO-AUTOInterface EMSO-AUTO

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20

1

2

3

4

5

6

7

8

9

10

x1

x 2

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.20

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x1

x 2

0.134 0.136 0.138 0.14 0.142 0.144 0.146 0.148

0.63

0.64

0.65

0.66

0.67

0.68

0.69

x1

x 2

Parameter p Eigenvalues Phase plane

p < 0.06361 Real negatives eigenvalues – stable nodep = 0.05 = [-1.13, -2.06]

p = 0.06361 Repeated real negatives eigenvalues – stable node (star) = [-1.4372, -1.4372]

0.06361 < p < 0.0889 Complex eigenvalues with negative real part - stable focusp = 0.085 = -1.095 0.565 i

36

Interface EMSO-AUTOInterface EMSO-AUTO

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20

1

2

3

4

5

6

7

8

9

10

x1

x 2

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

1

2

3

4

5

6

7

x1

x 2

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

1

2

3

4

5

6

7

x1

x 2

p = 0,0889unstable node gives rise to two points: unstable node and saddle point p > 0.0889

Turning point (fold): One stable solution (focus) and other unstable (node). (point 3 in figure below) = -1.009 0.605 i = [0, 3.432]

0.0889 < p < 0.0933 One stable solution (focus) and two unstable (saddle and node)p = 0.09 = -0.982 0.614 i = [-0.213, 3.332] = [0.364, 3.151]

0.0933 < p < 0.10574at p = 0.105738931 the first point goes from stable focus to stable node: = [-0.055, -0.046]

One stable solution (focus) and two unstable (saddle and focus)p = 0.10 = -0.652 0.651 i = [-0.439, 1.953] = 1.431 1.851 i

37

Interface EMSO-AUTOInterface EMSO-AUTO

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

1

2

3

4

5

6

7

x1

x 2

0 2 4 6 8 10 12 14 16 18 200

2

4

6

8

10

12

t

x

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20

2

4

6

8

10

12

x1

x 2

0.8 0.82 0.84 0.86 0.88 0.9 0.92 0.94 0.96 0.98 13

3.5

4

4.5

5

5.5

6

x1

x 2

p = 0.10574the stable node gives rise to two points: stable node and saddle forp < 0.10574

Turning point (fold): One stable solution (node), other unstable (focus), and one sable limit cycle. (point 2 in figure below) = [-0.097, 0] = 1.186 2.478 i

0.10574 < p < 0.1309 One unstable solution (focus) and one stable limit cyclep = 0.12 = 0.528 3.487 i

p = 0.1309 Hopf bifurcation: pure imaginary eigenvalues. (point 4 in figure below) = 4.008 i

38

Interface EMSO-AUTOInterface EMSO-AUTO

0.8 0.82 0.84 0.86 0.88 0.9 0.92 0.94 0.96 0.98 13

3.5

4

4.5

5

5.5

6

x1

x 2

p > 0.1309 Complex eigenvalues with negative real part - stable focusp = 0.15 = -0.952 4.627 i

p

x 2

39

Interface EMSO-AUTOInterface EMSO-AUTO

-4 -3 -2 -1 0 1 2 3 4-10

-8

-6

-4

-2

0

2

4

6

8

10

Re()

Im(

)

Hopf

Hopf

1st turning point

2nd turning point

Trajectories:stable pointsaddle pointunstable point

40

Interface EMSO-AUTOInterface EMSO-AUTO

Example: copy files auto_emso.exe and r-emso.bat (Windows) or @r-emso (linux) in “bin” folder of EMSO to the folder CSTR_auto and execute the command below in a prompt of commands (shell):

Windows: r-emso cstr_bif

Linux: ./@r-emso cstr_bif

The results are stored in file fort.7. In Linux the graphic tool PLAUT can be used to plot the results using the command @p.

41

Sensitivity AnalysisSensitivity Analysis

Objective: determine the effect of variation of parameters (p) or input variables (u) on the output variables.

Steady-state simulation: );,(

0);,(

puxHy

puxF

Sensitivity analysis

local:

global: bifurcation diagram, surface response

,

iy

j x p

yW

p

,

ix

j x p

xW

p

(case study)

1

x

F FW

x p

y x

H HW W

x p

,

j iy

i j x p

p yW

y p

Normalized form:

42

Sensitivity AnalysisSensitivity Analysis

Dynamic simulation: 0 0( , , , ; ) 0 , ( ; ) ( )

( , ; )

F t x x u p x t p x p

y H x u p

00( ) ( ) 0 , ( )x x x

xF F FW t W t W t

x x p p

p

HtW

x

HW xy

)(

where ( ) xx

dWW t

dt( )x

xW t

p

43

Sensitivity AnalysisSensitivity Analysis

44

Sensitivity AnalysisSensitivity Analysis

45

ReferencesReferences

DAE Solvers:

DASSL: Petzold, L.R. (1989), http://www.enq.ufrgs.br/enqlib/numeric/numeric.html

DASSLC: Secchi, A.R. (1992), http://www.enq.ufrgs.br/enqlib/numeric/numeric.html

MEBDFI: Abdulla, T.J. and J.R. Cash (1999), http://www.netlib.org/ode/mebdfi.f

PSIDE: Lioen, W.M., J.J.B. de Swart, and W.A. van der Veen (1997), http://www.cwi.nl/cwi/projects/PSIDE/

SUNDIALS: Serban, R. et al. (2004), http://www.llnl.gov/CASC/sundials/description/description.html