Operator representation of an networks of PDEs modeling a coronary - - PowerPoint PPT Presentation

operator representation of an networks of pdes modeling a
SMART_READER_LITE
LIVE PREVIEW

Operator representation of an networks of PDEs modeling a coronary - - PowerPoint PPT Presentation

Prirodoslovno-matemati cki fakultet Matemati cki odsjek Sveu cili ste u Zagrebu Operator representation of an networks of PDEs modeling a coronary stent OTIND 2016, Vienna, December 20, 2016 Luka Grubi si c, Josip Ivekovi


slide-1
SLIDE 1

Prirodoslovno-matematiˇ cki fakultet Matematiˇ cki odsjek Sveuˇ ciliˇ ste u Zagrebu

Operator representation of an networks of PDEs modeling a coronary stent

OTIND 2016, Vienna, December 20, 2016 Luka Grubiˇ si´ c, Josip Ivekovi´ c and Josip Tambaˇ ca Department of Mathematics, University of Zagreb luka.grubisic@math.hr

Luka Grubiˇ si´ c 1 / 44

slide-2
SLIDE 2

Outline

1 About stents

Description of a Stent

2 Mixed variational formulations 3 Stent modeling

1D model of curved rods 1D stent model

4 Eigenvalue problem for 1D stent model

Evolution problem Mixed formulation Finite element spaces on Graphs Examples

5 Some further developments

Luka Grubiˇ si´ c 2 / 44

slide-3
SLIDE 3

About stents

  • carotid stenosis

About stents Luka Grubiˇ si´ c 3 / 44

slide-4
SLIDE 4

About stents

  • carotid stenosis
  • stent: ”

solution”of the problem

About stents Luka Grubiˇ si´ c 3 / 44

slide-5
SLIDE 5

About stents Luka Grubiˇ si´ c 4 / 44

slide-6
SLIDE 6

About stents Luka Grubiˇ si´ c 4 / 44

slide-7
SLIDE 7

About stents Luka Grubiˇ si´ c 5 / 44

slide-8
SLIDE 8

Stent properties

  • A network of cylindrical tubes (made by laser cuts)
  • Typically: 316L stainless steel, but also cobalt, chrome and nickel.
  • expanded on the place of stenosis (balloon expandable is dominant

(99%) over self-expanding)

  • properties depend on

◮ complex geometry of a stent, ◮ mechanical properties of material.

  • metal

theory of elasticity

  • start from a conservation law

constrained evolution problems

About stents Luka Grubiˇ si´ c 6 / 44

slide-9
SLIDE 9

Plan

Mixed variational formulations are good for analyzing systems with constraints. Modal analysis gives first information on time evolution.

Mixed eigenvalue problems

Let us review two canonical approaches to mixed eigenvalue problems

About stents Luka Grubiˇ si´ c 7 / 44

slide-10
SLIDE 10

Type I

Dirichlet eigenvalue problem

−△u = λu u ∈ H = H1

  • No finite accumulation points
  • eigenfunctions make an orthonormal system for H.

A mixed formulation reads H = I − Div Grad

  • ,

M = −I

  • .

Modeli Luka Grubiˇ si´ c 8 / 44

slide-11
SLIDE 11

Type II

Stokes eigenvalue problem

Seek u ∈ H1

0([0, π]2) \ {0} and λ ∈ R such that

−△u + ∇p = λu ∇ · u = 0 .

  • no finite accumulation points
  • eigenfunctions span a subspace of H = L2(Ω)2.

A mixed formulation reads H = −△ Grad − Div

  • ,

M = I

  • .

Modeli Luka Grubiˇ si´ c 9 / 44

slide-12
SLIDE 12

Some pictures

IsoValue

  • 35.5992
  • 30.5136
  • 27.1232
  • 23.7328
  • 20.3424
  • 16.952
  • 13.5616
  • 10.1712
  • 6.7808
  • 3.3904
  • 1.68677e-008

3.3904 6.7808 10.1712 13.5616 16.952 20.3424 23.7328 27.1232 35.5992 Vec Value 0.0900351 0.18007 0.270105 0.360141 0.450176 0.540211 0.630246 0.720281 0.810316 0.900351 0.990386 1.08042 1.17046 1.26049 1.35053 1.44056 1.5306 1.62063 1.71067 Eigen Vector 0 valeur =52.3471

  • For Type I take P1–P0
  • For Typer II take P2–P1

Modeli Luka Grubiˇ si´ c 10 / 44

slide-13
SLIDE 13

Quasidefinite block operator matrices – form aproach

  • Split the sesquilinear form

h(u, v) = k(u1, v1) + b(u1, v2) + b(v1, u2), as h = k + b ≃ K B∗ B

  • .
  • representation theorems

self-adjoint operators

◮ Ker(h) = (Ker(K) ∩ Ker(B)) ⊕ Ker(B′) ◮ higher order Sobolev spaces and interpolation results? ◮ indefinite Cholesky on the level of block operator matrices e.g.

Grubiˇ si´ c, Kostrykin, Makarov, Veseli´ c 2013

  • This is indeed a special structure, recall the talk of C. Tretter or A.

Motovilov. With this we have access to spectral calculus to study constrained problems!

Modeli Luka Grubiˇ si´ c 11 / 44

slide-14
SLIDE 14

Representation theorems for quasidefinite matrices

For h we obtain a representation theorem h(u, v) = (Hu, v)H where H = K B∗ B

  • denotes the operator defined by h with

Dom(H) = (Dom(K) ∩ Dom(B)) ⊕ Dom(B∗) If K is invertible then H =

  • I

BK −1 I K −BK −1B∗ I K −1B∗ I

  • .

Has a flavor of a second representation theorem if the domain stability condition holds.

Modeli Luka Grubiˇ si´ c 12 / 44

slide-15
SLIDE 15

Similarity transformations

Since H =

  • I

BK −1 I K −BK −1B∗ I K −1B∗ I

  • .

and

  • I

BK −1 I −1 =

  • I

−BK −1 I

  • if follows

K −BK −1B∗ ˜ u ˜ v

  • = λ

−M ˜ u ˜ v

  • .

and so BK −1B∗˜ v = λM˜ v is an equivalent eigenvalue problem.

Modeli Luka Grubiˇ si´ c 13 / 44

slide-16
SLIDE 16

Similarity transformations

For Type II there are various equivalent reformulations (should K be singular) like K + B∗B B∗ B u p

  • = λ

M u p

  • .

since I B∗ I K B∗ B

  • =

K + B∗B B∗ B

  • and so M stays unchanged.

Recall H−1 = K −1 − K −1B∗(BK −1B∗)−1BK −1 K −1B∗(BK −1B∗)−1 (BK −1B∗)−1BK −1 −(BK −1B∗)−1

  • .

Modeli Luka Grubiˇ si´ c 14 / 44

slide-17
SLIDE 17

Finally ...

u p

  • = λ

(K −1 − K −1B∗(BK −1B∗)−1BK −1)Mu (BK −1B∗)−1BK −1Mu

  • .

and so u p

  • = λ

TMu SMu

  • .

which yields a reduced problem (in Ker(B)) T −1u = λMu . All finite eigenvalues of the reduced and the full problem coincide. Infinity is an eigenvalue of the saddle-point problem and it has associated vectors (cf. Spence et.al.).

Modeli Luka Grubiˇ si´ c 15 / 44

slide-18
SLIDE 18

Standard variational approach

Define solution operator T : H → Dom(k) and S : H → Dom(B), k(Tf , ˜ uS) + b(Sf , ˜ uS) = (f , ˜ uS), ˜ uS ∈ Dom(k), b(˜ nS, Tf ) = 0, ˜ nS ∈ Dom(B∗). (2.1) Require only ellipticity of kS in Ker(B), Recall Ker(H) = (Ker(K) ∩ Ker(B)) ⊕ Ker(B∗), eg. Benzi-Liesen-Golub, for matrices Tretter, Veseli´ c and Veseli´ c for operators. When K invertible T = (K −1 − K −1B∗(BK −1B∗)−1BK −1) further information cf. Boffi–Brezzi–Marini

Modeli Luka Grubiˇ si´ c 16 / 44

slide-19
SLIDE 19

Constraint satisfaction and computability

Direct computation yields BT = = B(K −1 − K −1B∗(BK −1B∗)−1BK −1) = BK −1 − BK −1B∗(BK −1B∗)−1BK −1 = (I − (BK −1B∗)(BK −1B∗)−1)BK −1 = 0 The action of the solution operator T is computable by solving the saddle point system.

Modeli Luka Grubiˇ si´ c 17 / 44

slide-20
SLIDE 20

Outline

1 About stents

Description of a Stent

2 Mixed variational formulations 3 Stent modeling

1D model of curved rods 1D stent model

4 Eigenvalue problem for 1D stent model

Evolution problem Mixed formulation Finite element spaces on Graphs Examples

5 Some further developments

Modeli Luka Grubiˇ si´ c 18 / 44

slide-21
SLIDE 21

Stent modeling

  • Consider the stent as a 3D elastic body

◮ geometry (net structure) ◮ elastic material: returns to the initial position after deformation λ, µ, E Stent modeling Luka Grubiˇ si´ c 19 / 44

slide-22
SLIDE 22

Stent modeling

  • Consider the stent as a 3D elastic body

◮ geometry (net structure) ◮ elastic material: returns to the initial position after deformation λ, µ, E ◮ for small deformation =

⇒ use linearized elasticity

◮ 3D stent simulation Stent modeling Luka Grubiˇ si´ c 19 / 44

slide-23
SLIDE 23

Stent modeling

  • Consider the stent as a 3D elastic body

◮ geometry (net structure) ◮ elastic material: returns to the initial position after deformation λ, µ, E ◮ for small deformation =

⇒ use linearized elasticity

◮ 3D stent simulation (very expensive!) Stent modeling Luka Grubiˇ si´ c 19 / 44

slide-24
SLIDE 24
  • a very complex 3D structure and computationally expensive

Stent modeling Luka Grubiˇ si´ c 20 / 44

slide-25
SLIDE 25
  • a very complex 3D structure and computationally expensive
  • stent struts are thin

Stent modeling Luka Grubiˇ si´ c 20 / 44

slide-26
SLIDE 26
  • a very complex 3D structure and computationally expensive
  • stent struts are thin
  • model struts with 1D model of curved rods

(rigorous justification of this step Jurak, Tambaˇ ca (1999), (2001))

  • model the whole stent as“a metric graph”of 1D rods!

(mathematical justification in nonlinear elasticity, Tambaˇ

ca, Velˇ ci´ c (2010), Griso (2010))

Stent modeling Luka Grubiˇ si´ c 20 / 44

slide-27
SLIDE 27

1D model of curved rods

˜ p′ + ˜ f = 0, ˜ q′ + t × ˜ p = 0, ˜ ω′ + QHQT˜ q = 0, ˜ ω(0) = ˜ ω(ℓ) = 0, ˜ u′ + t × ˜ ω = 0, ˜ u(0) = ˜ u(ℓ) = 0,

Stent modeling Luka Grubiˇ si´ c 21 / 44

slide-28
SLIDE 28

1D model of curved rods

˜ p′ + ˜ f = 0, ˜ q′ + t × ˜ p = 0, ˜ ω′ + QHQT˜ q = 0, ˜ ω(0) = ˜ ω(ℓ) = 0, ˜ u′ + t × ˜ ω = 0, ˜ u(0) = ˜ u(ℓ) = 0,

  • system od 12 ODE

Stent modeling Luka Grubiˇ si´ c 21 / 44

slide-29
SLIDE 29

1D model of curved rods

˜ p′ + ˜ f = 0, ˜ q′ + t × ˜ p = 0, ˜ ω′ + QHQT˜ q = 0, ˜ ω(0) = ˜ ω(ℓ) = 0, ˜ u′ + t × ˜ ω = 0, ˜ u(0) = ˜ u(ℓ) = 0,

  • system od 12 ODE
  • Q = (t, n, b) – rotation (Frenet basis for middle curve)
  • t – tangent, n – normal, b – binormal for middle curve

Stent modeling Luka Grubiˇ si´ c 21 / 44

slide-30
SLIDE 30

1D model of curved rods

˜ p′ + ˜ f = 0, ˜ q′ + t × ˜ p = 0, ˜ ω′ + QHQT˜ q = 0, ˜ ω(0) = ˜ ω(ℓ) = 0, ˜ u′ + t × ˜ ω = 0, ˜ u(0) = ˜ u(ℓ) = 0,

  • system od 12 ODE
  • Q = (t, n, b) – rotation (Frenet basis for middle curve)
  • t – tangent, n – normal, b – binormal for middle curve
  • H – matrix, properties of the cross-section of the rod, material

properties

Stent modeling Luka Grubiˇ si´ c 21 / 44

slide-31
SLIDE 31

1D model of curved rods

˜ p′ + ˜ f = 0, ˜ q′ + t × ˜ p = 0, ˜ ω′ + QHQT˜ q = 0, ˜ ω(0) = ˜ ω(ℓ) = 0, ˜ u′ + t × ˜ ω = 0, ˜ u(0) = ˜ u(ℓ) = 0,

  • system od 12 ODE
  • Q = (t, n, b) – rotation (Frenet basis for middle curve)
  • t – tangent, n – normal, b – binormal for middle curve
  • H – matrix, properties of the cross-section of the rod, material

properties

  • the equations are – balance of contact forces, balance of contact

couples, constitutive relationship for the material, in-extensibility and unsherability condition.

Stent modeling Luka Grubiˇ si´ c 21 / 44

slide-32
SLIDE 32

1D model of curved rods

˜ p′ + ˜ f = 0, ˜ q′ + t × ˜ p = 0, ˜ ω′ + QHQT˜ q = 0, ˜ ω(0) = ˜ ω(ℓ) = 0, ˜ u′ + t × ˜ ω = 0, ˜ u(0) = ˜ u(ℓ) = 0,

  • system od 12 ODE
  • Q = (t, n, b) – rotation (Frenet basis for middle curve)
  • t – tangent, n – normal, b – binormal for middle curve
  • H – matrix, properties of the cross-section of the rod, material

properties

  • the equations are – balance of contact forces, balance of contact

couples, constitutive relationship for the material, in-extensibility and unsherability condition.

  • 1D model much easier and quicker to solve than 3D

Stent modeling Luka Grubiˇ si´ c 21 / 44

slide-33
SLIDE 33

1D model of curved rods

˜ p′ + ˜ f = 0, ˜ q′ + t × ˜ p = 0, ˜ ω′ + QHQT˜ q = 0, ˜ ω(0) = ˜ ω(ℓ) = 0, ˜ u′ + t × ˜ ω = 0, ˜ u(0) = ˜ u(ℓ) = 0,

  • system od 12 ODE
  • Q = (t, n, b) – rotation (Frenet basis for middle curve)
  • t – tangent, n – normal, b – binormal for middle curve
  • H – matrix, properties of the cross-section of the rod, material

properties

  • the equations are – balance of contact forces, balance of contact

couples, constitutive relationship for the material, in-extensibility and unsherability condition.

  • 1D model much easier and quicker to solve than 3D

◮ numerical approximation in 3D: minutes ◮ numerical approximation in 1D: seconds Stent modeling Luka Grubiˇ si´ c 21 / 44

slide-34
SLIDE 34

1D stent model

  • At each edge: the system of 12 ODE

˜ pe′ + ˜ f

e = 0,

˜ qe′ + te × ˜ pe = 0, ˜ ωe′ + QeHe(Qe)T˜ qe = 0, ˜ ue′ + te × ˜ ωe = 0,

Stent modeling Luka Grubiˇ si´ c 22 / 44

slide-35
SLIDE 35

1D stent model

  • At each edge: the system of 12 ODE

˜ pe′ + ˜ f

e = 0,

˜ qe′ + te × ˜ pe = 0, ˜ ωe′ + QeHe(Qe)T˜ qe = 0, ˜ ue′ + te × ˜ ωe = 0,

  • Junction conditions

◮ for kinematical quantities: ω, u continuity Stent modeling Luka Grubiˇ si´ c 22 / 44

slide-36
SLIDE 36

1D stent model

  • At each edge: the system of 12 ODE

˜ pe′ + ˜ f

e = 0,

˜ qe′ + te × ˜ pe = 0, ˜ ωe′ + QeHe(Qe)T˜ qe = 0, ˜ ue′ + te × ˜ ωe = 0,

  • Junction conditions

◮ for kinematical quantities: ω, u continuity ◮ for dynamical quantities: p, q

1 sum of contact forces (˜

p) at each vertex =0

2 sum of contact couples (˜

q) at each vertex =0

Stent modeling Luka Grubiˇ si´ c 22 / 44

slide-37
SLIDE 37

1D stent model

  • At each edge: the system of 12 ODE

˜ pe′ + ˜ f

e = 0,

˜ qe′ + te × ˜ pe = 0, ˜ ωe′ + QeHe(Qe)T˜ qe = 0, ˜ ue′ + te × ˜ ωe = 0,

  • Junction conditions

◮ for kinematical quantities: ω, u continuity ◮ for dynamical quantities: p, q

1 sum of contact forces (˜

p) at each vertex =0

2 sum of contact couples (˜

q) at each vertex =0

rigorous mathematical justification in nonlinear elasticity

(Tambaˇ ca, Velˇ ci´ c (2010), Griso (2010))

Stent modeling Luka Grubiˇ si´ c 22 / 44

slide-38
SLIDE 38

1D stent model

Unknown: U = (U1, . . . , UnE) = ((˜ u1, ˜ ω1), . . . , (˜ unE, ˜ ωnE)) Test function: V = (V1, . . . , VnE) = ((˜ v1, ˜ w1), . . . , (˜ vnE, ˜ wnE))

Stent modeling Luka Grubiˇ si´ c 23 / 44

slide-39
SLIDE 39

1D stent model

Unknown: U = (U1, . . . , UnE) = ((˜ u1, ˜ ω1), . . . , (˜ unE, ˜ ωnE)) Test function: V = (V1, . . . , VnE) = ((˜ v1, ˜ w1), . . . , (˜ vnE, ˜ wnE)) Function spaces: H1(N; R6) =

  • V ∈

nE

  • i=1

H1((0, ℓe); R6) : Vi((Φi)−1(v)) = Vj((Φj)−1(v)), v ∈ V, v ∈ ei ∩ ej , Vstent ={V ∈ H1(N; R6) : ˜ vi ′ + ti × ˜ wi = 0, i = 1, . . . , nE,

  • N

V = 0}.

Stent modeling Luka Grubiˇ si´ c 23 / 44

slide-40
SLIDE 40

1D stent model

Unknown: U = (U1, . . . , UnE) = ((˜ u1, ˜ ω1), . . . , (˜ unE, ˜ ωnE)) Test function: V = (V1, . . . , VnE) = ((˜ v1, ˜ w1), . . . , (˜ vnE, ˜ wnE)) Function spaces: H1(N; R6) =

  • V ∈

nE

  • i=1

H1((0, ℓe); R6) : Vi((Φi)−1(v)) = Vj((Φj)−1(v)), v ∈ V, v ∈ ei ∩ ej , Vstent ={V ∈ H1(N; R6) : ˜ vi ′ + ti × ˜ wi = 0, i = 1, . . . , nE,

  • N

V = 0}. Sum up weak formulations for rods:

Stent modeling Luka Grubiˇ si´ c 23 / 44

slide-41
SLIDE 41

1D stent model

Unknown: U = (U1, . . . , UnE) = ((˜ u1, ˜ ω1), . . . , (˜ unE, ˜ ωnE)) Test function: V = (V1, . . . , VnE) = ((˜ v1, ˜ w1), . . . , (˜ vnE, ˜ wnE)) Function spaces: H1(N; R6) =

  • V ∈

nE

  • i=1

H1((0, ℓe); R6) : Vi((Φi)−1(v)) = Vj((Φj)−1(v)), v ∈ V, v ∈ ei ∩ ej , Vstent ={V ∈ H1(N; R6) : ˜ vi ′ + ti × ˜ wi = 0, i = 1, . . . , nE,

  • N

V = 0}. Sum up weak formulations for rods: find U ∈ Vstent such that

nE

  • i=1

ℓi QiHi(Qi)T(˜ ωi)′ · ˜ w′dx1 =

nE

  • i=1

ℓi ˜ f

i · ˜

vdx1, V ∈ Vstent

(Tambaˇ ca, Kosor, ˇ Cani´ c, Paniagua, SIAM J. Appl. Math., 2010)

Stent modeling Luka Grubiˇ si´ c 23 / 44

slide-42
SLIDE 42

1D stent model: weak formulation revisited

a(U, V) =

nE

  • i=1

ℓi QiHi(Qi)T(˜ ωi)′ · ˜ w′dx1, b(V, Ξ) =

nE

  • i=1

ℓi ξi · (˜ vi ′ + ti × ˜ wi)dx1 + α ·

nE

  • i=1

ℓi ˜ vidx1 + β ·

nE

  • i=1

ℓi ˜ widx1, Ξ = (ξ1, . . . , ξnE, α, β), f (V) =

nE

  • i=1

ℓi ˜ f

i · ˜

vidx1, M = L2(N; R3) × R3 × R3 =

nE

  • i=1

L2(0, ℓi; R3) × R3 × R3 K = Vstent = {V ∈ H1(N; R6) : b(V, Θ) = 0, ∀Θ ∈ M}.

Stent modeling Luka Grubiˇ si´ c 24 / 44

slide-43
SLIDE 43

1D stent model: weak formulation revisited

a(U, V) =

nE

  • i=1

ℓi QiHi(Qi)T(˜ ωi)′ · ˜ w′dx1, b(V, Ξ) = unsherability/inextensibility +the nonslip condition. Ξ = (ξ1, . . . , ξnE, α, β), f (V) =

nE

  • i=1

ℓi ˜ f

i · ˜

vidx1, M = L2(N; R3) × R3 × R3 =

nE

  • i=1

L2(0, ℓi; R3) × R3 × R3 K = Vstent = {V ∈ H1(N; R6) : b(V, Θ) = 0, ∀Θ ∈ M}.

Stent modeling Luka Grubiˇ si´ c 25 / 44

slide-44
SLIDE 44

Evolution problem

˜ pi ′ + ˜ f

i = ρiAi∂tt˜

ui, ˜ qi ′ + ti × ˜ pi = 0, ˜ ωi ′ + QiHi(Qi)T˜ qi = 0, ˜ ui ′ + ti × ˜ ωi = θi, i = 1, . . . , nE + junction conditions

Eigenvalue problem for 1D stent model Luka Grubiˇ si´ c 26 / 44

slide-45
SLIDE 45

Evolution problem

˜ pi ′ + ˜ f

i = ρiAi∂tt˜

ui, ˜ qi ′ + ti × ˜ pi = 0, ˜ ωi ′ + QiHi(Qi)T˜ qi = 0, ˜ ui ′ + ti × ˜ ωi = θi, i = 1, . . . , nE + junction conditions Let (limit of 3D linearized Antman-Cosserat model) m(U, V) =

nE

  • i=1

ρiAi ℓi ˜ ui · ˜ vidx1,

Problem (EvoP)

Find U ∈ L2(0, T; Vstent) such that d2 dt2 m(U, V) + a(U, V) = f (V), V ∈ Vstent.

Eigenvalue problem for 1D stent model Luka Grubiˇ si´ c 26 / 44

slide-46
SLIDE 46

Convergence theory by Boffi–Brezzi–Marini

Recall solution operator T : H → Dom(k) and S : H → Dom(B), k(Tf , ˜ uS) + b(Sf , ˜ uS) = (f , ˜ uS), ˜ uS ∈ Dom(k), b(˜ nS, Tf ) = 0, ˜ nS ∈ Dom(B∗). (4.2)

Disscretization

For a discretization take Xh i Mh as piecewise polynomial spaces on N. Then the discrete solution operators are k(Thfh, ˜ uS) + b(Shf , ˜ uS) = (f , ˜ uS), ˜ uS ∈ Xh, b(˜ nS, Tf ) = 0, ˜ nS ∈ Mh. (4.3) The space Kh ⊂ K is the space ov polynomials which satisfy the constraints.

Eigenvalue problem for 1D stent model Luka Grubiˇ si´ c 27 / 44

slide-47
SLIDE 47

Convergence theory by Boffi–Brezzi–Gastaldi

Proposition

Let Kh be a nonempty sequence of subspaces for which we have inperpolation estimates and let k be elliptic on K. Then there is ω(h) = o(1) such that Tf − Thf VS ≤ ω(h)f L2(N;Z3). Bounded compact operators Th norm converge to T. If the resolvent is converging somewhere – say at z = 0 –then it converges for every z in resolvent set ρ(T) = C \ Spec(T).

Eigenvalue problem for 1D stent model Luka Grubiˇ si´ c 28 / 44

slide-48
SLIDE 48

Convergence rate – coments

Let λ, u and λh and uh be eigenvalues and eigenvectors from Xh Then |λ − λh| ≤ CT − ThDom(k) = O(ω(h)) u − uhDom(k) ≤ CT − ThDom(k) = O(ω(h)) . eigenfuction c.rate optimal, eigenvalue c.rate not. If Sh exists, then |λ − λh| ≤ CT − ThDom(k)S − Sh . λh is a Ritz value of h but not a Ritz value of the solution operator T.

Eigenvalue problem for 1D stent model Luka Grubiˇ si´ c 29 / 44

slide-49
SLIDE 49

Finite element spaces

0.005 0.01 0.015 0.02 −2 −1 1 2 x 10

−3

−1.5 −1 −0.5 0.5 1 1.5 x 10

−3

  • For finite element spaces take

X (n)

h

= {uS ∈ H1(N; R6) : uS|ei ∈ Pn(ei; R6), i = 1, . . . , nE}, M(m)

h

= {nS ∈ L2(N; R6) : nS|ei ∈ Pm(ei; R3), i = 1, . . . , nE}.

  • We expect the error behaving as in 1D interpolation.
  • Piecewise polynomial approximation of the curved middle line?

Eigenvalue problem for 1D stent model Luka Grubiˇ si´ c 30 / 44

slide-50
SLIDE 50

Convergence rates

  • Sobolev spaces on graphs – nice review by O. Post
  • Good interpolation operators for lower order spaces – hard because of

the contact conditions in junctions! Geometry of the graph plays a role.

  • For X (n)

h , M(n−1) h

approximation we obtain (#DOF)−(n+1) decay rates for the eigenvalue errors.

  • for second order problems see Arioli and Benzi 2015.

The interplay of geometry and constraints

Doing interpolation on each edge and then assembling into a graph fails for first order polynomials since the constraint is to restrictive. Using higher order polynomials introduces numerical integration problems in matrix assembly.

Eigenvalue problem for 1D stent model Luka Grubiˇ si´ c 31 / 44

slide-51
SLIDE 51

Matrix formulation

For stiffness we have K =       · · · A11 A12 · · · A1m+1 . . . . . . . . . . . . ... . . . . . . · · · Am+11 Am+12 · · · Am+1m+1       , B =        B11 · · · B1n+1 C 11 · · · C 1n+1 . . . · · · . . . Bm+11 · · · Bm+1n+1 C m+11 · · · C m+1n+1        . and for mass M =

  • M
  • =

        D11 D12 · · · D1m+1 · · · . . . . . . . . . . . . ... . . . . . . Dm+11 Dm+12 · · · Dm+1m+1 · · ·         .

  • translations have mass but no stifness
  • micro-rotations have stifness but no mass
  • graph connectivity is in B
  • stifness positive definite on K.

Eigenvalue problem for 1D stent model Luka Grubiˇ si´ c 32 / 44

slide-52
SLIDE 52

Examples of eigenproblem for 1D stent model

Four stents considered

0.005 0.01 0.015 0.02 −2 −1 1 2 x 10

−3

−1.5 −1 −0.5 0.5 1 1.5 x 10

−3

0.01 0.02 −2 −1 1 2 x 10

−3

−1.5 −1 −0.5 0.5 1 1.5 x 10

−3

0.005 0.01 0.015 0.02 −2 −1 1 2 x 10

−3

−1.5 −1 −0.5 0.5 1 1.5 x 10

−3

0.01 0.02 −1.5 −1 −0.5 0.5 1 1.5 x 10

−3

−1.5 −1 −0.5 0.5 1 1.5 x 10

−3

Eigenvalue problem for 1D stent model Luka Grubiˇ si´ c 33 / 44

slide-53
SLIDE 53

Leading eigenvalues

Palmaz Cypher Express Xience 1. 1.033 0.8894 0.06014 0.05488 2. 1.033 0.8895 0.06014 0.05488 3. 5.265 1.3683 0.32504 0.28767 4. 7.499 3.5328 0.33972 0.32201 5. 7.499 3.5329 0.33973 0.32201 6. 11.329 3.6604 0.58740 0.58038

Eigenvalue problem for 1D stent model Luka Grubiˇ si´ c 34 / 44

slide-54
SLIDE 54

Palmaz

−5 5 10 15 20 x 10

−3

−3 −2 −1 1 2 3 x 10

−3

−5 5 10 15 20 x 10

−3

−3 −2 −1 1 2 3 x 10

−3

−5 5 10 15 20 x 10

−3

−1.5 −1 −0.5 0.5 1 1.5 x 10

−3

−5 5 10 15 20 x 10

−3

−3 −2 −1 1 2 3 x 10

−3

−5 5 10 15 20 x 10

−3

−3 −2 −1 1 2 3 x 10

−3

−2 −1.5 −1 −0.5 0.5 1 1.5 2 x 10

−3

−2 −1 1 2 x 10

−3

Eigenvalue problem for 1D stent model Luka Grubiˇ si´ c 35 / 44

slide-55
SLIDE 55

Cypher

−5 5 10 15 20 x 10

−3

−3 −2 −1 1 2 3 x 10

−3

−5 5 10 15 20 x 10

−3

−3 −2 −1 1 2 3 x 10

−3

0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 −2 −1 1 2 x 10

−3

−5 5 10 15 20 x 10

−3

−3 −2 −1 1 2 3 x 10

−3

−5 5 10 15 20 x 10

−3

−3 −2 −1 1 2 3 x 10

−3

−4 −3 −2 −1 1 2 3 4 x 10

−3

−1.5 −1 −0.5 0.5 1 1.5 x 10

−3

Eigenvalue problem for 1D stent model Luka Grubiˇ si´ c 36 / 44

slide-56
SLIDE 56

Express

−5 5 10 15 20 x 10

−3

−3 −2 −1 1 2 3 x 10

−3

−5 5 10 15 20 x 10

−3

−3 −2 −1 1 2 3 x 10

−3

−5 5 10 15 20 x 10

−3

−2 −1 1 2 x 10

−3

−5 5 10 15 20 x 10

−3

−3 −2 −1 1 2 3 x 10

−3

−5 5 10 15 20 x 10

−3

−3 −2 −1 1 2 3 x 10

−3

−5 5 10 15 20 x 10

−3

−2 −1 1 2 x 10

−3

Eigenvalue problem for 1D stent model Luka Grubiˇ si´ c 37 / 44

slide-57
SLIDE 57

Xience

−5 5 10 15 20 x 10

−3

−3 −2 −1 1 2 3 x 10

−3

−5 5 10 15 20 x 10

−3

−3 −2 −1 1 2 3 x 10

−3

−5 5 10 15 20 x 10

−3

−2 −1 1 2 x 10

−3

−5 5 10 15 20 x 10

−3

−3 −2 −1 1 2 3 x 10

−3

−5 5 10 15 20 x 10

−3

−3 −2 −1 1 2 3 x 10

−3

−5 5 10 15 20 x 10

−3

−2 −1 1 2 x 10

−3

Eigenvalue problem for 1D stent model Luka Grubiˇ si´ c 38 / 44

slide-58
SLIDE 58

Convergence rates — dominated by geometry errors

Eigenvalue problem for 1D stent model Luka Grubiˇ si´ c 39 / 44

slide-59
SLIDE 59

Exponential integrators

Recall T −1u = ω2Mu . so let’s use, for u0 ∈ Ker(B), u = cos((TM)−1/2 t)u0 + (TM)1/2 sin((TM)−1/2 t)u1 for vibration analysis.

Some further developments Luka Grubiˇ si´ c 40 / 44

slide-60
SLIDE 60

Exponential integrators

Recall T −1u = ω2Mu . so let’s use, for u0 ∈ Ker(B), u = cos((TM)−1/2 t)u0 + (TM)1/2 sin((TM)−1/2 t)u1 for vibration analysis. For the“heat”equation see Emmrich and Mehrmann

Some further developments Luka Grubiˇ si´ c 40 / 44

slide-61
SLIDE 61

Exponential integrators

Recall T −1u = ω2Mu . so let’s use, for u0 ∈ Ker(B), u = cos((TM)−1/2 t)u0 + (TM)1/2 sin((TM)−1/2 t)u1 for vibration analysis. For the“heat”equation see Emmrich and Mehrmann Grimm, Hochbruck, Goeckler (trigonometric integrators) use resolvent Krylov subspace projection to propagate a very stiff second order system.

Some further developments Luka Grubiˇ si´ c 40 / 44

slide-62
SLIDE 62

Trigonometric integrators

Action of T is accessible by solving the saddle-point problem. For u0, so that u0 + Gu0 < ∞ we have u − Wn(τG 1/2)u0 ≤ τ C √n Gu0 where Wn is the rational approximation of the sine or cosine operator from the (resolvent) Krylov space Kn = span{Tu0, T 2u0, · · · , T nu0} 1D stent simulation

Some further developments Luka Grubiˇ si´ c 41 / 44

slide-63
SLIDE 63

Some conclusions

  • Can detect undesirable vibration modes (meshy structures are stiffer,

but can unexpectedly buckle).

  • For a numerical integration we would like to approximate curved rod

by splitting it in a chain of straight rods which approximate the middle line

  • Questions ?

◮ Discrete inf-sup ◮ Good linear for Schur complements? Some further developments Luka Grubiˇ si´ c 42 / 44

slide-64
SLIDE 64

Acknowledgment

Research supported by the Croatian Science Foundation grant nr. HRZZ 9345

Some further developments Luka Grubiˇ si´ c 43 / 44

slide-65
SLIDE 65

Thank you for your attention!

Some further developments Luka Grubiˇ si´ c 44 / 44