An application of numerical bifurcation analysis Greg Lewis - - PowerPoint PPT Presentation

An application of numerical bifurcation analysis

Greg Lewis

University of Ontario Institute of Technology (UOIT)

with Bill Langford (Guelph) and Wayne Nagata (UBC)

BIRS, August 8, 2007

BIRS – p.1/48

Outline

Introduction The differentially heated rotating annulus experiment Bifurcation analysis Numerical continuation Eigenvalue computation Examples differentially heated rotating annulus differentially heated rotating spherical shell Summary

BIRS – p.2/48

A differentially heated rotating annulus

D fluid r

b

T

a

T

b a

r Ω R

y r x φ z

Ω

BIRS – p.3/48

A differentially heated rotating planet

South Pole North Pole

cold hot

BIRS – p.4/48

A differentially heated rotating annulus

D fluid r

b

T

a

T

b a

r Ω R

y r x φ z

Ω

BIRS – p.5/48

Regime diagram

lower symmetric vacillation irregular steady waves upper symmetric knee

log(Thermal Rossby number) log(Taylor number)

BIRS – p.6/48

Wave flow in the annulus

BIRS – p.7/48

Vacillating flow in the annulus

BIRS – p.8/48

Regime diagram

lower symmetric vacillation irregular steady waves upper symmetric knee

log(Thermal Rossby number) log(Taylor number)

BIRS – p.9/48

Bifurcation analysis

Nonlinear DE: dx

dt = G(x, α), x ∈ Rn, α ∈ R1.

Steady solution x0 = x0(α) when: G(x0, α) = 0. Look for bifurcations from steady solution linear stability of steady solution from eigenvalues, λ, of the linearization of dynamical equation about the steady solution:

Gx(x = x0, α). Real (λj) < 0 for all j → x0 is linearly stable Real (λj) > 0 for one j → x0 is linearly unstable

BIRS – p.10/48

Numerical computations

Steady solutions use pseudo-arclength continuation Linear stability: eigenvalues Implicitly restarted Arnoldi method with Cayley transformations

BIRS – p.11/48

Steady solution: continuation

Look for steady solutions discretization reduces PDE to system of nonlinear algebraic equations need to solve G(x, α) = 0, x ∈ Rn, α ∈ R Use Newton’s method with continuation need to have a good guess assume we know x0 at α0 such that G(x0, α0) = 0

BIRS – p.12/48

Natural parameterization x α

G(x, )=0 α

α0

α1

x x

1

x1

^

BIRS – p.13/48

Natural parameterization

x0 x1

α G(x, ) = 0

x1

α 0 α

x

α 1

BIRS – p.14/48

Pseudo-arclength continuation

Consider the parameter α as an unknown predictor: new guess (ˆ

x1, ˆ α1) given by ˆ x1 = x0 + ∆s t0t(x)

0 ,

ˆ α1 = α0 + ∆s t0t(α) t0 = [t(x) t(α)

0 ] is the tangent to the solution curve

the step size ∆s measures arclength along tangent line for corrector, add an extra condition to get new system:

G(x, α) = f(x, α) =

BIRS – p.15/48

Pseudo-arclength continuation

x0 x1

α G(x, ) = 0

x2 x3

t 2 t 1

x2 x3

α 0 α

x

α 1

BIRS – p.16/48

Eigenvalue approximation

Eigenvalue problem Linearize about steady solution get generalized eigenvalue problems

λBΦ = AΦ

discretization leads to matrix eigenvalue problems

BIRS – p.17/48

Eigenvalue approximation

For eigenvalues use ‘Implicitly restarted Arnoldi method’ iterative memory efficient finds extremal eigenvalues

BIRS – p.18/48

Eigenvalue approximation

Use generalized Cayley transform

C(A, B) = (A − σ1B)−1 (A − σ2B) λ are eigenvalues from λBx = Ax µ are eigenvalues from µx′ = Cx′ Real(λ) > σ1 + σ2 2 → |µ| > 1

BIRS – p.19/48

Eigenvalue approximation

Use generalized Cayley transform

C(A, B) = (A − σ1B)−1 (A − σ2B)

Don’t need to form the matrix C explicitly

• nly need the matrix-vector product w = Cv

w = Cv = (A − σ1B)−1 (A − σ2B) v

multiple by (A − σ1B) get:

(A − σ1B) w = (A − σ2B) v

i.e. a system of linear equations

BIRS – p.20/48

Centre manifold reduction

Apply centre manifold reduction at bifurcation points gives a low-dimensional model of dynamics get existence and stability of bifurcating solutions gives results close to a bifurcation point (local dynamics) Write ODE (reduced equation) in normal form compute the coefficients of the normal form equations Deduce dynamics of PDE from low-dimensional ODE

BIRS – p.21/48

A differentially heated rotating annulus

D fluid r

b

T

a

T

b a

r Ω R

y r x φ z

Ω

BIRS – p.22/48

Model of fluid in the annulus

Navier-Stokes equations in the Boussinesq approximation Cylindrical coordinates and rotating frame of reference No-slip boundary conditions Insulating top and bottom of annulus Differential heating: ∆T = Tb − Ta inner cylinder cooled; outer cylinder heated Quantitatively accurate results

BIRS – p.23/48

Analysis

Look for steady flows invariant under rotation primary transitions reduces to problem in two-spatial dimensions Bifurcations from steady solutions

BIRS – p.24/48

Regime diagram

lower symmetric vacillation irregular steady waves upper symmetric knee

log(Thermal Rossby number) log(Taylor number)

BIRS – p.25/48

Transition curve

10

5

10

6

10

7

10

8

10

−3

10

−2

10

−1

10

Taylor number thermal Rossby number

(3,4) (4,5) (5,6) (6,7) (7,8) (8,7) (7,6) (6,5) theoretical transition curve theoretical critical wave number transitions experimental transition curve experimental critical wave number transitions

BIRS – p.26/48

Regions of bi-stability

10

5

10

6

10

7

10

8

10

−3

10

−2

10

−1

10 Taylor number thermal Rossby number (3,4) (4,5) (5,6) (6,7) (7,8) (8,7) (7,6) (6,5) theoretical transition curve theoretical critical wave number transitions boundaries of region of bistability

BIRS – p.27/48

Spherical Shell

∆T

g Ω

BIRS – p.28/48

Model of fluid in a spherical shell

Navier-Stokes equations in the Boussinesq approximation Spherical polar coordinates and rotating frame of reference No-slip boundary conditions at inner sphere Stress-free boundary condition at outer sphere Insulating outer sphere Differential heating imposed on inner sphere: at r = r0, T = T0 − ∆T cos(2θ).

BIRS – p.29/48

Differential heating

0.5 1 1.5 1.5 2 2.5 3 3.5 4 4.5

θ T T = T0 T = T0 − ∆ T cos(2 θ) pole equator

BIRS – p.30/48

Spherical shell

r0 R

∆T

Ω

BIRS – p.31/48

Analysis

Look for steady flows invariant under rotation and reflection about equator Reduces to problem in two-spatial dimensions Introduces additional boundary conditions at pole and equator Bifurcations of steady solutions

BIRS – p.32/48

Steady Solution:

,

x y stream function 1 2 0.5 1 1.5 2 x y azimuthal fluid velocity

+ −

1 2 0.5 1 1.5 2 x y temperature deviation

− +

1 2 0.5 1 1.5 2

Bristol 05 – p.10/24

Steady Solution:

,

x y stream function 1 2 0.5 1 1.5 2 x y azimuthal fluid velocity

+ −

1 2 0.5 1 1.5 2 x y temperature deviation

− +

1 2 0.5 1 1.5 2

Bristol 05 – p.11/24

Steady Solution:

,

x y stream function 1 2 0.5 1 1.5 2 x y azimuthal fluid velocity

+ −

1 2 0.5 1 1.5 2 x y temperature deviation

− +

1 2 0.5 1 1.5 2

Bristol 05 – p.12/24

Steady Solution:

,

x y stream function 1 2 0.5 1 1.5 2 x y azimuthal fluid velocity

+ −

1 2 0.5 1 1.5 2 x y temperature deviation

− +

1 2 0.5 1 1.5 2

Bristol 05 – p.13/24

Steady Solution:

,

x y stream function 1 2 0.5 1 1.5 2 x y azimuthal fluid velocity

+ −

1 2 0.5 1 1.5 2 x y temperature deviation

− +

1 2 0.5 1 1.5 2

Bristol 05 – p.14/24

Bifurcation Diagram:

0.01 0.0125 0.015 0.0175 0.02 −9.7 −9.6 x 10

−5

Eigenvalue with largest real part ∆ T max real(λ)

0.01 0.0125 0.015 0.0175 0.02 3 4 5 6 x 10

−3

Continuation of steady solution 1−cell 2−cell ∆ T || ξ ||2

Bristol 05 – p.15/24

Steady Solution:

,

x y stream function 1 2 0.5 1 1.5 2 x y azimuthal fluid velocity

+ −

1 2 0.5 1 1.5 2 x y temperature deviation

− +

1 2 0.5 1 1.5 2

Bristol 05 – p.16/24

Steady Solution:

,

x y stream function 1 2 0.5 1 1.5 2 x y azimuthal fluid velocity

+ −

1 2 0.5 1 1.5 2 x y temperature deviation

− +

1 2 0.5 1 1.5 2

Bristol 05 – p.17/24

Bifurcation Diagram:

0.0107 0.0108 0.0109 −10 −5 x 10

−8

Eigenvalue with largest real part ∆ T max real(λ)

0.0107 0.0108 0.0109 1 1.5 x 10

−3

Continuation of steady solution 1−cell 2−cell ∆ T || ξ ||2

DEDS: Pattern Formation – p.27/35

Cusp bifurcation

∆T ∆T ∆T

small η η = η0 large η

HP SN SN

Dynamics Day – p.28/35

Cusp bifurcation (schematic)

∆T

3 solutions 1 solution

η

BIRS – p.43/48

Computation of cusp point

Codimension two bifurcation Need two parameters: ∆T and η Write equations as:

˙ U = LU + N(U, U)

where U is dependent variable,

LU is linear part, N(U, U) is nonlinear part,

and ˙

U is derivative with respect to time

BIRS – p.44/48

Computation of cusp point

Cusp point is characterized by:

• 1. LU0 + N(U0, U0) = 0
• 2. zero eigenvalue of L0 where

L0V = LV + N(V, U0) + N(U0, V )

• 3. vanishing of the coefficient of 2nd-order term of

equation on centre manifold (or reduced equation)

BIRS – p.45/48

Reduced equation

Reduced equation

˙ w = β1 + β2w + aw2 + cw3

where

a = 1/2 Φ∗, N(Φ, Φ) = 0 Φ is the eigenfunction corresponding to λ = 0, Φ∗ is the corresponding adjoint eigenfunction, ·, · is the inner product

BIRS – p.46/48

Defining system

LU0 + N(U0, U0) = 0, g = 0, g′ = 0

where g and g′ are scalars given by

L0V + gB = 0, C, V = 1 L0V ′ + g′B = −N(V, V ),

• C, V ′

= 0

where B not in range of L0, and C not in range of the adjoint operator L∗

0.

Solve to get a = 0 at η = 3.46, ∆T = 0.011

BIRS – p.47/48

Summary

Application of numerical bifurcation analysis compute flow regimes compute details of flow transitions Could apply same ideas to industrial problems Applied to transitions from steady flows Could also apply similar ideas to transitions from periodic flows HPC

BIRS – p.48/48