DIANA An object-oriented tool for nonlinear analysis of chemical - - PowerPoint PPT Presentation

2 September 2008

MAX−PLANCK−INSTITUT TECHNISCHER SYSTEME MAGDEBURG DYNAMIK KOMPLEXER

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

DIANA — An object-oriented tool for nonlinear analysis

• f chemical processes

Mykhaylo Krasnyk

Max Planck Institute for Dynamics of Complex Technical Systems, PSPD group Otto-von-Guericke-University, IFAT

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Outline

1

Introduction

2

Steady-state point analysis Limit points continuation Higher co-dimension singularities Simulation models in Diana Parameter continuation

3

Case Study I: Nonlinear analysis of MCFC

4

Periodic solutions continuation Analysis of periodic solutions Recursive Projection Method

5

Case Study II: Periodic solutions in MSMPR Crystallizer

6

Summary

OvGU, IFAT Outline 2/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Introduction

Chemical Process Engineering

OvGU, IFAT Introduction 3/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Introduction

Chemical Process Engineering Dynamical Systems Analysis

cin ˙ q c T, c

OvGU, IFAT Introduction 3/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Introduction

Chemical Process Engineering Dynamical Systems Analysis

cin ˙ q c T, c t [s] T [K]

OvGU, IFAT Introduction 3/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Introduction

Chemical Process Engineering Dynamical Systems Analysis

cin ˙ q c T, c t [s] T [K] ˙ q [l/h] T [K]

OvGU, IFAT Introduction 3/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Introduction

Chemical Process Engineering Dynamical Systems Analysis

cin ˙ q c T, c t [s] T [K] ˙ q [l/h] T [K]

State of the art

AUTO/MATCONT detection and continuation of bifurcation points and pe- riodic solutions in low-order ODE, no modeling tool DIVA stationary and dynamic simulations of higher-order DAE for engineering processes, outdated FORTRAN77 code LOCA bifurcation analysis of large-scale CFD applications, lim- ited application problems

OvGU, IFAT Introduction 3/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Software tool Diana

Diana — Dynamic simulation and nonlinear analysis tool

developed at MPI Magdeburg modularization, extensibility and object-oriented architecture equation based models numerical solvers based on free code enhanced scripting and visualization

OvGU, IFAT Introduction 4/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Software tool Diana

Diana — Dynamic simulation and nonlinear analysis tool

developed at MPI Magdeburg modularization, extensibility and object-oriented architecture equation based models numerical solvers based on free code enhanced scripting and visualization

Objectives of the work

generation of C++ model code for Diana symbolic differentiation of models (Maxima package) parameter continuation of nonlinear problems

higher-order singularities of steady-state curves efficient calculation of periodic solutions by reduction techniques

analysis of test models

OvGU, IFAT Introduction 4/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Steady-state point continuation

Parameter continuation vs. dynamic simulation

Condition for steady-state points of an autonomous system x ∈ Rn, ν ∈ Rp is: ˙ x = f (x, ν)

!

=0

OvGU, IFAT Steady-state point analysis 5/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Steady-state point continuation

Parameter continuation vs. dynamic simulation

Condition for steady-state points of an autonomous system x ∈ Rn, ν ∈ Rp is: ˙ x = f (x, ν)

!

=0

x t xs x(t, x0)

OvGU, IFAT Steady-state point analysis 5/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Steady-state point continuation

Parameter continuation vs. dynamic simulation

Condition for steady-state points of an autonomous system x ∈ Rn, ν ∈ Rp is: ˙ x = f (x, ν)

!

=0

x t xs x(t, x0)

→

x λ ∈ ν x(t, x0)|t=∞ xs λs

OvGU, IFAT Steady-state point analysis 5/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Steady-state point continuation

Parameter continuation vs. dynamic simulation

Condition for steady-state points of an autonomous system x ∈ Rn, ν ∈ Rp is: ˙ x = f (x, ν)

!

=0

x t xs x(t, x0)

→

x λ ∈ ν x(λ) xs λs

OvGU, IFAT Steady-state point analysis 5/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Steady-state point continuation

Parameter continuation vs. dynamic simulation

Condition for steady-state points of an autonomous system x ∈ Rn, ν ∈ Rp is: ˙ x = f (x, ν)

!

=0

x t xs x(t, x0)

→

x λ ∈ ν x(λ) xs λs

Stability is determined by eigenvalues of a linearized system at the steady-state point xs, νs: V Λ = ∂ f ∂ x V

OvGU, IFAT Steady-state point analysis 5/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Limit points analysis

Example: limit and hysteresis points (α, β, λ ∈ ν)

x λ

❙ ▲ ▲ ❍ ▲ ▲ ❍ P

OvGU, IFAT Steady-state point analysis/ Limit points continuation 6/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Limit points analysis

Example: limit and hysteresis points (α, β, λ ∈ ν)

x λ

❙ ▲ ▲ ❍ ▲ ▲ ❍ P

OvGU, IFAT Steady-state point analysis/ Limit points continuation 6/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Limit points analysis

Example: limit and hysteresis points (α, β, λ ∈ ν)

x λ α

❙ ▲ ▲ ❍ ▲ ▲ ❍ P

OvGU, IFAT Steady-state point analysis/ Limit points continuation 6/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Limit points analysis

Example: limit and hysteresis points (α, β, λ ∈ ν)

x λ α

❙ ▲ ▲ ❍

projection to α-λ plane

α λ

▲ ▲ ❍ P

OvGU, IFAT Steady-state point analysis/ Limit points continuation 6/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Limit points analysis

Example: limit and hysteresis points (α, β, λ ∈ ν)

x λ α

❙ ▲ ▲ ❍

projection to α-λ plane

α λ β

▲ ▲ ❍ P

OvGU, IFAT Steady-state point analysis/ Limit points continuation 6/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Lyapunov-Schmidt reduction 1

Reduction definition

For a system ˙ x = f (x, ν) with f (xs, νs) = 0 and L = fx(xs, νs) with dim ker L = 1 analysis of a limit points curve can be performed with a scalar equation g(z, ν)!

• 1M. Golubitsky and D. G. Schaeffer, Singularities and groups in bifurcation theory. Vol. I, 1985.

OvGU, IFAT Steady-state point analysis/ Limit points continuation 7/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Lyapunov-Schmidt reduction 1

Reduction definition

For a system ˙ x = f (x, ν) with f (xs, νs) = 0 and L = fx(xs, νs) with dim ker L = 1 analysis of a limit points curve can be performed with a scalar equation g(z, ν)! Lyapunov-Schmidt reduction defines a mapping φ : ker L → range L⊥ φ(v, ν) := (I − E)f (v + W (v, ν), ν)

• 1M. Golubitsky and D. G. Schaeffer, Singularities and groups in bifurcation theory. Vol. I, 1985.

OvGU, IFAT Steady-state point analysis/ Limit points continuation 7/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Lyapunov-Schmidt reduction 1

Reduction definition

For a system ˙ x = f (x, ν) with f (xs, νs) = 0 and L = fx(xs, νs) with dim ker L = 1 analysis of a limit points curve can be performed with a scalar equation g(z, ν)! Lyapunov-Schmidt reduction defines a mapping φ : ker L → range L⊥ φ(v, ν) := (I − E)f (v + W (v, ν), ν) The basis v0 ∈ ker L and v ∗

0 ∈ range L⊥ is defined by an adjoint system

8 > < > : f (x, ν) = 0, fx(x, ν)v0 − βv ∗ = 0, ||v0||2 = 1, f T

x (x, ν)v ∗ 0 − γv0 = 0,

||v ∗

0 ||2 = 1.

Reduced equation in the chosen basis {v0, v ∗

0 } is

g(z, λ) = v ∗

0 , f (zv0 + W (zv0, λ), λ),

where z ∈ R and λ ∈ R ⊂ ν

• 1M. Golubitsky and D. G. Schaeffer, Singularities and groups in bifurcation theory. Vol. I, 1985.

OvGU, IFAT Steady-state point analysis/ Limit points continuation 7/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Higher co-dimension singularities

Classification of singular points with codim g 3

The classification theorem guarantees existence only the following possible singu- larities of g with codim g 3

codim g 1 2 3 z5 ± λ z4 ± zλ z3 ± λ2 z2 ± λ4 z4 ± λ

gzzzz = 0 gλ = 0

z3 ± zλ

gzzz = 0 gzλ = 0

z2 ± λ3

gzz = 0 | d3g| = 0

z3 ± λ

gzzz = 0 gλ = 0

z2 ± λ2

gzz = 0 | d2g| = 0

z2 ± λ

gzz = 0 gλ = 0

equilibrium

gz = 0

type

OvGU, IFAT Steady-state point analysis/ Higher co-dimension singularities 8/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Parameter continuation

Parameter continuation solver

Predictor-corrector method computes a parametrized curve y(ζ) := {x(ζ), ν(ζ)} such that F(y(ζ)) ≡ 0 with 0 ζ ζmax — arc-length parameter

OvGU, IFAT Steady-state point analysis/ Parameter continuation 9/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Parameter continuation

Parameter continuation solver

Predictor-corrector method computes a parametrized curve y(ζ) := {x(ζ), ν(ζ)} such that F(y(ζ)) ≡ 0 with 0 ζ ζmax — arc-length parameter In Diana two types of predictors are used: tangential predictor y (k)

t x ν y(k)

• y (k)

t

y(ζ)

OvGU, IFAT Steady-state point analysis/ Parameter continuation 9/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Parameter continuation

Parameter continuation solver

Predictor-corrector method computes a parametrized curve y(ζ) := {x(ζ), ν(ζ)} such that F(y(ζ)) ≡ 0 with 0 ζ ζmax — arc-length parameter In Diana two types of predictors are used: tangential predictor y (k)

t

chord predictor y (k)

c x ν y(k−1) y(k)

• y (k)

c

y(ζ)

OvGU, IFAT Steady-state point analysis/ Parameter continuation 9/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Parameter continuation

Parameter continuation solver

Predictor-corrector method computes a parametrized curve y(ζ) := {x(ζ), ν(ζ)} such that F(y(ζ)) ≡ 0 with 0 ζ ζmax — arc-length parameter Types of correctors are used: local parametrization y (k+1)

i

= ˜ y (k+1)

i x ν y(k) ˜ y(k+1) y(k+1) y(ζ)

OvGU, IFAT Steady-state point analysis/ Parameter continuation 9/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Parameter continuation

Parameter continuation solver

Predictor-corrector method computes a parametrized curve y(ζ) := {x(ζ), ν(ζ)} such that F(y(ζ)) ≡ 0 with 0 ζ ζmax — arc-length parameter Types of correctors are used: local parametrization y (k+1)

i

= ˜ y (k+1)

i

pseudo-arclength parametrization y (k+1) − ˜ y (k+1) ⊥ y (k+1)

x ν y(k) ˜ y(k+1) y(k+1) y(ζ)

OvGU, IFAT Steady-state point analysis/ Parameter continuation 9/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Outline

1

Introduction

2

Steady-state point analysis Limit points continuation Higher co-dimension singularities Simulation models in Diana Parameter continuation

3

Case Study I: Nonlinear analysis of MCFC

4

Periodic solutions continuation Analysis of periodic solutions Recursive Projection Method

5

Case Study II: Periodic solutions in MSMPR Crystallizer

6

Summary

OvGU, IFAT Case Study I: Nonlinear analysis of MCFC 10/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Model: Molten Carbonate Fuel Cell

Model equations

H Ly Lz ΦAE ΦA ΦCE ΦC anode electrolyte cathode Itot ∆Φtot

Energy balance and the corresponding bound- ary conditions are ∂Θ ∂τ = ∂2Θ ∂η2 + ` B − φtot´ i − Bi1Θ, ∂Θ ∂η ˛ ˛ ˛ ˛

0,τ

= Bi2Θ(0, τ), ∂Θ ∂η ˛ ˛ ˛ ˛

1,τ

= −Bi2Θ(1, τ) A correlation for i due to the Butler-Volmer reaction kinetics is i = ψA exp „ γA Θ 1 + Θ « ( exp −(1 − βA)γeq φA 1 + Θ ! − K A

eq exp

βAγeq φA 1 + Θ !) Ohm’s law for the electrolyte is i = ψE exp „ γE Θ 1 + Θ « “ φA + φC − φtot” , Itot =

1

Z i dη The simulation model is discretized with an equidistant grid and with 100 grid points and has 201 variables

OvGU, IFAT Case Study I: Nonlinear analysis of MCFC 11/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Results of the MCFC analysis

Unfolding results near a winged cusp point

1 1.5 2 2.5 3 10 12 14 16

Bi2 γA Pitchfork continuation

OvGU, IFAT Case Study I: Nonlinear analysis of MCFC 12/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Results of the MCFC analysis

Unfolding results near a winged cusp point

1 1.5 2 2.5 3 10 12 14 16

Bi2 γA Pitchfork continuation

z5 ± λ z4 ± zλ z3 ± λ2 z2 ± λ4 z4 ± λ z3 ± xλ z2 ± λ3 z3 ± λ z2 ± λ2 z2 ± λ equilibrium

OvGU, IFAT Case Study I: Nonlinear analysis of MCFC 12/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Results of the MCFC analysis

Unfolding results near a winged cusp point

1 1.5 2 2.5 3 10 12 14 16

Bi2 γA

• Pitchfork continuation

10 10.5 11 11.5 12 12.5 0.95 1 1.05

γA γE

• codim-1 singularities

z5 ± λ z4 ± zλ z3 ± λ2 z2 ± λ4 z4 ± λ z3 ± xλ z2 ± λ3 z3 ± λ z2 ± λ2 z2 ± λ equilibrium

OvGU, IFAT Case Study I: Nonlinear analysis of MCFC 12/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Results of the MCFC analysis

Unfolding results near a winged cusp point

1 1.5 2 2.5 3 10 12 14 16

Bi2 γA

• Pitchfork continuation

10 10.5 11 11.5 12 12.5 0.95 1 1.05

γA γE

• codim-1 singularities

z5 ± λ z4 ± zλ z3 ± λ2 z2 ± λ4 z4 ± λ z3 ± xλ z2 ± λ3 z3 ± λ z2 ± λ2 z2 ± λ equilibrium

OvGU, IFAT Case Study I: Nonlinear analysis of MCFC 12/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Results of the MCFC analysis

Unfolding results near a winged cusp point

1 1.5 2 2.5 3 10 12 14 16

Bi2 γA

• Pitchfork continuation

10 10.5 11 11.5 12 12.5 0.95 1 1.05

γA γE

• codim-1 singularities

z5 ± λ z4 ± zλ z3 ± λ2 z2 ± λ4 z4 ± λ z3 ± xλ z2 ± λ3 z3 ± λ z2 ± λ2 z2 ± λ equilibrium

OvGU, IFAT Case Study I: Nonlinear analysis of MCFC 12/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Results of the MCFC analysis

Unfolding results near a winged cusp point

1 1.5 2 2.5 3 10 12 14 16

Bi2 γA

• Pitchfork continuation

10 10.5 11 11.5 12 12.5 0.95 1 1.05

γA γE

• codim-1 singularities

1 2 3 4 5 ×105 10 20 30 40

Itot Θ0 Steady-state point continuation

z5 ± λ z4 ± zλ z3 ± λ2 z2 ± λ4 z4 ± λ z3 ± xλ z2 ± λ3 z3 ± λ z2 ± λ2 z2 ± λ equilibrium

OvGU, IFAT Case Study I: Nonlinear analysis of MCFC 12/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Results of the MCFC analysis

Unfolding results near a winged cusp point

1 1.5 2 2.5 3 10 12 14 16

Bi2 γA

• Pitchfork continuation

10 10.5 11 11.5 12 12.5 0.95 1 1.05

γA γE

• codim-1 singularities

1 2 3 4 ×104 1 2

Itot Θ0 Steady-state point continuation

z5 ± λ z4 ± zλ z3 ± λ2 z2 ± λ4 z4 ± λ z3 ± xλ z2 ± λ3 z3 ± λ z2 ± λ2 z2 ± λ equilibrium

OvGU, IFAT Case Study I: Nonlinear analysis of MCFC 12/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Results of the MCFC analysis

Unfolding results near a winged cusp point

1 1.5 2 2.5 3 10 12 14 16

Bi2 γA

• Pitchfork continuation

10 10.5 11 11.5 12 12.5 0.95 1 1.05

γA γE

• codim-1 singularities

1 2 3 4 ×104 1 2

Itot Θ0 Steady-state point continuation

z5 ± λ z4 ± zλ z3 ± λ2 z2 ± λ4 z4 ± λ z3 ± xλ z2 ± λ3 z3 ± λ z2 ± λ2 z2 ± λ equilibrium

OvGU, IFAT Case Study I: Nonlinear analysis of MCFC 12/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Results of the MCFC analysis

Unfolding results near a winged cusp point

1 1.5 2 2.5 3 10 12 14 16

Bi2 γA

• Pitchfork continuation

10 10.5 11 11.5 12 12.5 0.95 1 1.05

γA γE

• codim-1 singularities

1 2 3 4 ×104 1 2

Itot Θ0 Steady-state point continuation

0.25 0.5 0.75 1 5 10 15

η Θ(η) Spatial profile, Itot = 37000

OvGU, IFAT Case Study I: Nonlinear analysis of MCFC 12/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Outline

1

Introduction

2

Steady-state point analysis Limit points continuation Higher co-dimension singularities Simulation models in Diana Parameter continuation

3

Case Study I: Nonlinear analysis of MCFC

4

Periodic solutions continuation Analysis of periodic solutions Recursive Projection Method

5

Case Study II: Periodic solutions in MSMPR Crystallizer

6

Summary

OvGU, IFAT Periodic solutions continuation 13/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Task statement

Fixed point problem formulation

Find a fixed point {x∗ ∈ Rn, γ∗ ∈ Rm} of a mapping F with a constraint G x(i+1) = F(x(i), γ), F : Rn × Rm → Rn ∈ C ∞, n 1, 0 = G(x, γ), G : Rn × Rm → Rm ∈ C ∞, such that x∗ = F(x∗, γ∗) and G(x∗, γ∗) = 0.

OvGU, IFAT Periodic solutions continuation/ Analysis of periodic solutions 14/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Task statement

Fixed point problem formulation

Find a fixed point {x∗ ∈ Rn, γ∗ ∈ Rm} of a mapping F with a constraint G x(i+1) = F(x(i), γ), F : Rn × Rm → Rn ∈ C ∞, n 1, 0 = G(x, γ), G : Rn × Rm → Rm ∈ C ∞, such that x∗ = F(x∗, γ∗) and G(x∗, γ∗) = 0.

Calculation of periodic solutions

For an IVP with ϕ(t, x0, λ), such that f (ϕ, ˙ ϕ, λ) ≡ 0 F(x, γ) := ϕ(x, T, λ) — Poincar´ e map, G(x, γ) :=  s(x, T, λ) n(x, T, λ) ff x(0) = x0 ∈ Rn γ = {T, λ} ∈ R2 with the pseudo-arclength parameterisation n(x, γ) := x − ˜ x, ˜ xt + γ − ˜ γ, ˜ γt − σ and the phase condition s(x, T, λ) := Z T ϕ(t, x, λ), ˙ ˜ ϕ(t) dt

t ˜ ϕ ϕ

OvGU, IFAT Periodic solutions continuation/ Analysis of periodic solutions 14/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

RPM: Idea of the method2

The stabilization procedure

Let M := F ∗

x have eigenvalues {µk}n 1 that for some δ > 0 are ordered as

|µ1| · · · |µnp| > 1 − δ |µnp+1| · · · |µn| and np n

• 2G. M. Shroff and H. B. Keller. SIAM JNA, 30(4):1099–1120, Aug. 1993.

OvGU, IFAT Periodic solutions continuation/ Recursive Projection Method 15/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

RPM: Idea of the method2

The stabilization procedure

Let M := F ∗

x have eigenvalues {µk}n 1 that for some δ > 0 are ordered as

|µ1| · · · |µnp| > 1 − δ |µnp+1| · · · |µn| and np n Then possible to define the maximal invariant subspace U of M belonging to {µk}

np 1

with projectors P and Q := I − P that induce an orthogonal direct sum decompo- sition Rn = U ⊕ U ⊥ = PRn ⊕ QRn

• 2G. M. Shroff and H. B. Keller. SIAM JNA, 30(4):1099–1120, Aug. 1993.

OvGU, IFAT Periodic solutions continuation/ Recursive Projection Method 15/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

RPM: Idea of the method2

The stabilization procedure

Let M := F ∗

x have eigenvalues {µk}n 1 that for some δ > 0 are ordered as

|µ1| · · · |µnp| > 1 − δ |µnp+1| · · · |µn| and np n Then possible to define the maximal invariant subspace U of M belonging to {µk}

np 1

with projectors P and Q := I − P that induce an orthogonal direct sum decompo- sition Rn = U ⊕ U ⊥ = PRn ⊕ QRn A subspace decomposition of the original system leads to 2 6 4 V T

q F (i) x Vq − Iq

V T

p F (i) x Vp − Ip V T p F (i) γ

G (i)

x Vp

G (i)

γ

3 7 5 2 4 ∆¯ q(i) ∆¯ p(i) ∆γ(i) 3 5 = − 2 6 4 V T

q (r (i) + F (i) γ ∆γ(i))

V T

p (r (i) + F (i) x Vq∆¯

q(i)) G (i) + G (i)

x Vq∆¯

q(i) 3 7 5

• 2G. M. Shroff and H. B. Keller. SIAM JNA, 30(4):1099–1120, Aug. 1993.

OvGU, IFAT Periodic solutions continuation/ Recursive Projection Method 15/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

RPM: Idea of the method2

The stabilization procedure

Let M := F ∗

x have eigenvalues {µk}n 1 that for some δ > 0 are ordered as

|µ1| · · · |µnp| > 1 − δ |µnp+1| · · · |µn| and np n Then possible to define the maximal invariant subspace U of M belonging to {µk}

np 1

with projectors P and Q := I − P that induce an orthogonal direct sum decompo- sition Rn = U ⊕ U ⊥ = PRn ⊕ QRn A subspace decomposition of the original system leads to 2 6 4 V T

q F (i) x Vq − Iq

V T

p F (i) x Vp − Ip V T p F (i) γ

G (i)

x Vp

G (i)

γ

3 7 5 2 4 ∆¯ q(i) ∆¯ p(i) ∆γ(i) 3 5 = − 2 6 4 V T

q (r (i) + F (i) γ ∆γ(i))

V T

p (r (i) + F (i) x Vq∆¯

q(i)) G (i) + G (i)

x Vq∆¯

q(i) 3 7 5 The resulting system has properties: eigenvalues of the restricted to U⊥ matrix are |λ(V T

q F (i) x Vq)| 1 − δ

relatively small size of the bottom right part (np + m) × (np + m)

• 2G. M. Shroff and H. B. Keller. SIAM JNA, 30(4):1099–1120, Aug. 1993.

OvGU, IFAT Periodic solutions continuation/ Recursive Projection Method 15/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

RPM: Algorithm

Start P subspace calculation Q part solution P part solution Convergence? End no yes

Initial data

The input values are: starting values x(0) and γ(0) basis for the P part V (0)

p

= [v1, . . . , vnp+npe] Power subspace iteration approximates dominant np + npe eigenvalues of F (i)

x

(Arnoldi iteration) Find ∆q with an iterative method (Picard, Jacobi, Gauss-Seidel iteration, GMRES) Solve with a direct method one Newton step for ∆¯ p and ∆γ Results of the RPM algorithm: the values x∗ and γ∗ np dominant eigenvalues µ and corresponding eigenvectors Vp of the mapping derivative F ∗

x

OvGU, IFAT Periodic solutions continuation/ Recursive Projection Method 16/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

RPM: Algorithm

Start P subspace calculation Q part solution P part solution Convergence? End no yes

P subspace calculation

The input values are: starting values x(0) and γ(0) basis for the P part V (0)

p

= [v1, . . . , vnp+npe] Power subspace iteration approximates dominant np + npe eigenvalues of F (i)

x

(Arnoldi iteration) Find ∆q with an iterative method (Picard, Jacobi, Gauss-Seidel iteration, GMRES) Solve with a direct method one Newton step for ∆¯ p and ∆γ Results of the RPM algorithm: the values x∗ and γ∗ np dominant eigenvalues µ and corresponding eigenvectors Vp of the mapping derivative F ∗

x

OvGU, IFAT Periodic solutions continuation/ Recursive Projection Method 16/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

RPM: Algorithm

Start P subspace calculation Q part solution P part solution Convergence? End no yes

Q part solution (Picard iteration, ith step)

The input values are: starting values x(0) and γ(0) basis for the P part V (0)

p

= [v1, . . . , vnp+npe] Power subspace iteration approximates dominant np + npe eigenvalues of F (i)

x

(Arnoldi iteration) Find ∆q with an iterative method (Picard, Jacobi, Gauss-Seidel iteration, GMRES) Solve with a direct method one Newton step for ∆¯ p and ∆γ Results of the RPM algorithm: the values x∗ and γ∗ np dominant eigenvalues µ and corresponding eigenvectors Vp of the mapping derivative F ∗

x

OvGU, IFAT Periodic solutions continuation/ Recursive Projection Method 16/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

RPM: Algorithm

Start P subspace calculation Q part solution P part solution Convergence? End no yes

P part solution (Newton iteration, ith step)

The input values are: starting values x(0) and γ(0) basis for the P part V (0)

p

= [v1, . . . , vnp+npe] Power subspace iteration approximates dominant np + npe eigenvalues of F (i)

x

(Arnoldi iteration) Find ∆q with an iterative method (Picard, Jacobi, Gauss-Seidel iteration, GMRES) Solve with a direct method one Newton step for ∆¯ p and ∆γ Results of the RPM algorithm: the values x∗ and γ∗ np dominant eigenvalues µ and corresponding eigenvectors Vp of the mapping derivative F ∗

x

OvGU, IFAT Periodic solutions continuation/ Recursive Projection Method 16/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

RPM: Algorithm

Start P subspace calculation Q part solution P part solution Convergence? End no yes

Convergence check and results

The input values are: starting values x(0) and γ(0) basis for the P part V (0)

p

= [v1, . . . , vnp+npe] Power subspace iteration approximates dominant np + npe eigenvalues of F (i)

x

(Arnoldi iteration) Find ∆q with an iterative method (Picard, Jacobi, Gauss-Seidel iteration, GMRES) Solve with a direct method one Newton step for ∆¯ p and ∆γ Results of the RPM algorithm: the values x∗ and γ∗ np dominant eigenvalues µ and corresponding eigenvectors Vp of the mapping derivative F ∗

x

OvGU, IFAT Periodic solutions continuation/ Recursive Projection Method 16/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Outline

1

Introduction

2

Steady-state point analysis Limit points continuation Higher co-dimension singularities Simulation models in Diana Parameter continuation

3

Case Study I: Nonlinear analysis of MCFC

4

Periodic solutions continuation Analysis of periodic solutions Recursive Projection Method

5

Case Study II: Periodic solutions in MSMPR Crystallizer

6

Summary

OvGU, IFAT Case Study II: Periodic solutions in MSMPR Crystallizer 17/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Model: MSMPR Crystallizer 3

Model equations

Feed q, cin qq, c, F q, c, Fhp Product classification qf , c, Fhf Fines dissolver

The population balance equation is ∂F ∂t = −∂(GF) ∂L − q V (hf (L) + hp(L))F with the boundary condition F(0, t) = B(c, t) G(c, t) = kb(c(t) − csat)b kg(c(t) − csat)g . The recycle ratio of the fines dissolution and the classified product removal are hf = R1(1 − h(L − Lf )), hp = 1 + R2h(L − Lp) The mass balance of solute is MA dc dt = (ρ − MAc) „ q V + 1 ε dε dt « + qMAcin V ε − qρ V ε „ 1 + kv Z ∞ (hp − 1)FL3 dL « , where ε is the void fraction which is given by ε = 1 − kv Z ∞ FL3 dL

• 3P. K. Pathath and A. Kienle. CES,57:4391–4399,Oct. 2002.

OvGU, IFAT Case Study II: Periodic solutions in MSMPR Crystallizer 18/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Continuation results

Results of analysis (R2 = 0)

Hopf point continuation

5 10 15 20 50 100

b [-] R1

II I

OvGU, IFAT Case Study II: Periodic solutions in MSMPR Crystallizer 19/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Continuation results

Results of analysis (R2 = 0)

Hopf point continuation

5 10 15 20 50 100

b [-] R1

5 10 15 20 4.1 4.2 4.3 4.4

b [-] c [mol/l]

R1 = 0

OvGU, IFAT Case Study II: Periodic solutions in MSMPR Crystallizer 19/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Continuation results

Results of analysis (R2 = 0)

Hopf point continuation

5 10 15 20 50 100

b [-] R1

5 10 15 20 4.1 4.2 4.3 4.4

b [-] c [mol/l]

R1 = 0

5 10 15 20 4.1 4.2 4.3 4.4

b [-] c [mol/l]

R1 = 75

OvGU, IFAT Case Study II: Periodic solutions in MSMPR Crystallizer 19/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Continuation results

Results of analysis (R2 = 0)

Hopf point continuation

5 10 15 20 50 100

b [-] R1

5 10 15 20 4.1 4.2 4.3 4.4

b [-] c [mol/l]

R1 = 0

5 10 15 20 4.1 4.2 4.3 4.4

b [-] c [mol/l]

R1 = 75

5 10 15 20 4.1 4.2 4.3 4.4

b [-] c [mol/l]

R1 = 94

OvGU, IFAT Case Study II: Periodic solutions in MSMPR Crystallizer 19/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Continuation results

Results of analysis (R2 = 0)

Hopf point continuation

5 10 15 20 50 100

b [-] R1

5 10 15 20 4.1 4.2 4.3 4.4

b [-] c [mol/l]

R1 = 0

5 10 15 20 4.1 4.2 4.3 4.4

b [-] c [mol/l]

R1 = 75

5 10 15 20 4.1 4.2 4.3 4.4

b [-] c [mol/l]

R1 = 94

5 10 15 20 4.1 4.2 4.3 4.4

b [-] c [mol/l]

R1 = 98

OvGU, IFAT Case Study II: Periodic solutions in MSMPR Crystallizer 19/20

O O T T V O N G U E R I C K E U N I V E R S I T Ä T M A G D E B U R G

Summary

Topic: Simulation system Diana

Lisp-module for the modeling tool ProMoT has been developed (transformation and symbolic differentiation of ProMoT models, C++ code generator) C++ and Python interaction with CAPE-OPEN interfaces solver classes for linear and nonlinear problems have been implemented

Topic: Numerical nonlinear analysis

implementation of continuation methods for steady-state, limit and Hopf points generation of adjoint systems and reduced test functions for singular points periodic solution continuation with simple bifurcations detection recursive projection method has been applied to speed-up computation

Topic: Model analysis

more detailed analysis of a molten carbonate fuel cell has been performed results of periodic solutions continuation of an MSMPR crystallizer has been confirmed with reasonable time speed-up

OvGU, IFAT Summary 20/20