Modeling and discretization methods for the numerical simulation of elastic stents [original]

Modeling and discretization methods for the numerical simulation

of elastic stents ∗

Luka Grubiˇsi´c†Matko Ljulj†Volker Mehrmann‡Josip Tambaˇca†

January 22, 2019

Abstract

A new model description for the numerical simulation of elastic stents is proposed. Based

on the new formulation an inf-sup inequality for the finite element discretization is proved

and the proof of the inf-sup inequality for the continuous problem is simplified. The new

formulation also leads to faster simulation times despite an increased number of variables.

The techniques also simplify the analysis and numerical solution of the evolution problem

describing the movement of the stent under external forces. The results are illustrated via

numerical examples.

Keywords: elastic stent, mathematical modeling, numerical simulation, mixed finite element

formulation, stationary system, evolution equation,

AMS: 74S05, 74K10, 74K30, 74G15, 74H15, 65M15, 65M60

1 Introduction

In this paper we present a new model description for the dynamic and stationary simulation

of stents. This new formulation is using constrained partial differential equations in mixed

variational weak form, which are based on a network structure consisting of one dimensional

curved rods (struts). The new formulation will turn out to be particularly convenient for the

analysis of the partial differential equation, in particular in proving an inf-sup inequality which

directly transfers to a discrete inf-sup inequality in the discretized setting, so that from classical

results of [3] the error estimates follow.

In the new formulation, the inextensibility and unshearability of the rod are expressed in

the weak formulation, but the continuity of the displacement and the modeling of infinitesimal

rotations are modelled via constraints so that they do not have to be incorporated in the function

spaces as was done in the classical approach in [10]. This advantage of the new formulation comes

along with the introduction of new unknowns for the displacements and infinitesimal rotations

at vertices where different struts are connected and further unknowns for the contact couples

and contact forces at the end points of each strut. However, despite the introduction of many

new unknowns, the numerical solvers become more efficient for the large scale cases.

Stents, see Figure 1, are typically considered as a union of struts each of which is modeled by

a 1D curved rod model, see [14, 15], and a set of junction conditions describing the connection

∗This work has been supported by Deutscher Akademischer Austauschdienst (DAAD) via Project Asymptotic

and algebraic analysis of nonlinear eigenvalue problems in contact mechanics and electro magnetism. The third

author is also supported by Einstein Foundation Berlin via Einstein Center ECMath Project:Model Reduction

for Nonlinear Parameter-Dependent Eigenvalue Problems in Photonic Crystals.

†Department of Mathematics, Faculty of Science, University of Zagreb, Bijeniˇcka 30, 10000 Zagreb,

{luka,mljulj,tambaca}@math.hr.

§Institut f¨ur Mathematik MA 4-5, TU Berlin, Str. des 17. Juni 136, D-10623 Berlin, FRG.

[email protected].

Figure 1: Example of an elastic stent (Cypher stent by Cordis Corporation)

of the struts, see [9]. The so-obtained model describes the three-dimensional behavior of stents

but it has the complexity of a one-dimensional model. This model can be applied for any

elastic structure made of thin curved (or straight) rods. It was first formulated in [23] and then

reformulated in the weak form in [5]. The properties of the mixed formulation for the model

have been analyzed in [10], numerical methods have been introduced, and error estimates have

been derived in [12], but up to now error estimates for the contact forces, which are represented

as Lagrange multipliers of the contact conditions, were missing. To derive these error estimates

is one of the main results of this paper.

To model the topology of the stent we use an undirected graph N= (V,E) consisting of a

set Vof nVvertices, which are the points where the middle lines of the rods meet and a set Eof

nEedges that represent a 1D description of the curved rod. To be able to use a 1D curved rod

model, we additionally need to prescribe the local geometry of the rod, i.e., the middle curve

and the geometry of the cross-section as well as the material properties of the stent. These are

given by

•the function Φi: [0, `i]→R3as natural parametrization of the middle line of the ith strut

of length `i, represented by the edge ei∈ E,

•the shear modulus µiand the Young modulus Eias parameters describing the material of

the ith strut,

•as well as the width wiand the thickness tiof the rectangular cross-section of the ith strut.

Using these quantities, in the stationary case, see [23], the model for the ith strut ei∈ E, is

given by the following system of ordinary differential equations (in space)

0 = ∂spi+fi,(1.1)

0 = ∂sqi+ti×pi,(1.2)

0 = ∂sωi−Qi(Hi)−1(Qi)Tqi,(1.3)

0 = ∂sui+ti×ωi.(1.4)

where for the ith strut

•ui: [0, `i]→R3denotes the vector of displacements on the middle curve,

•ωi: [0, `i]→R3is the vector of infinitesimal rotations of the cross-section,

•qiis the contact moment and piis the contact force,

•fiis the line density of the applied forces,

•Qi= [ti,ni,bi] is an orthogonal rotation matrix associated to the middle curve, with

ti= (Φi)0being the unit tangent to the middle curve and ni,bibeing vectors spanning

the normal plane to the middle curve, so that Qirepresents the local basis at each point

of the middle curve,

•Hi= diag(µiKi, EiIi

n, EiIi

b) is a positive definite diagonal matrix, with the Young modulus

Ei, the shear modulus µi,Ii

n,Ii

bare the moments of inertia of the cross section and µiKi

is the torsional rigidity of the cross section.

Equations (1.1) and (1.2) represent equilibrium equations (for forces and moments), while (1.3)

and (1.4) are constitutive relations. In particular, (1.4) describes the inextensibility and un-

shearability of the struts, see [5] for more details.

In addition to equations (1.1)–(1.4), at each vertex of the stent we have a kinematic coupling

condition that uand ωare continuous and a dynamic coupling condition describing the balance

of contact forces pand contact moments q.

Denoting by J−

jthe set of all edges that leave the jth vertex, i.e., the local variable is equal

to 0 at vertex jand by J+

jthe set of all edges that enter the vertex, i.e., the local variable is

equal to `ifor ith edge at the vertex j. With these notations we obtain the node conditions

ωi(0) = ωk(`k), i ∈J−

j, k ∈J+

j, j = 1, . . . , nV,

ui(0) = uk(`k), i ∈J−

j, k ∈J+

j, j = 1, . . . , nV,

i∈J+

pi(`i)−X

i∈J−

pi(0) = 0, j = 1, . . . , nV,

i∈J+

qi(`i)−X

i∈J−

qi(0) = 0, j = 1, . . . , nV.

(1.5)

Since this is a pure traction problem, we can integrate over s∈[0, `i] and specify a unique

solution by requiring the two additional conditions

u:=

i=1 Z`i

uids = 0,ZN

ω:=

i=1 Z`i

ωids = 0,(1.6)

which means that the total displacement as well as the total infinitesimal rotation are zero.

In the formulation of [10], the model is described on the collection of all displacements ui

and infinitesimal rotations ωifor all edges which are continuous on the whole stent. Thus, the

tuples of unknowns in the problem uS= ((u1,ω1),...,(unE,ωnE)) belong to the space

H1(N;Rk) = (y1,...,ynE)∈

i=1

H1(0, `i;Rk) :

yi(0) = yk(`k), i ∈J−

j, k ∈J+

j, j = 1, . . . , nV,

with H1(0, `i;Rk) being the Sobolov space of functions on [0, `i] whose derivatives up to the

first derivative are square Lebesgue integrable. The formulation given in [10] is a mixed for-

mulation with Lagrange multipliers appearing in the formulation due to the inextensibility and

unshearability of the struts in the 1d curved rod model (1.4), and the two conditions on the

total displacement and infinitesimal rotation (1.6). Under these conditions, in [10] an inf-sup

inequality was proved and the well-posedness of the problem was established. However, for the

discretized problem via the finite element method it would be necessary to also have a discrete

inf-sup inequality to obtain an error estimate also for the Lagrange multipliers approximation.

We will show that the elastic energy is coercive in the space implementing inextensibility,

unshearability, and continuity of displacement and infinitesimal rotation. Using the results of [3],

we will show that the Lagrange multipliers are unique in both the old and the new formulation

and for the discrete problem also in the new formulation. This is a considerable improvement

compared to the classical mixed formulation [12], where a proof of the discrete inf-sup inequality

was not successful.

Let AI∈R3nV,3nEdenote the incidence matrix of the oriented graph (V,E) with three

connected components, organized in the following way: a 3×3 submatrix at rows 3i−2,3i−1,3i

and columns 3j−2,3j−1,3jis I3if the edge jenters the vertex i,−I3if it leaves the vertex

ior 0 otherwise. Then the matrix A+

Iis obtained from AIby setting all elements −1 to 0 and

A−

Iis defined as AI=A+

I−A−

I. Let us also introduce the projectors

E∈R3,3nE,Pj

V∈R3,3nV

on the coordinates 3i−2,3i−1,3iand 3j−2,3j−1,3j, respectively. We will also need the

spaces

L2(N;R3) =

i=1

L2(0, `i;R3), L2

Hr(N;R3) =

i=1

Hr(0, `i;R3), r ≥1

with associated norms

k(y1,...,ynE)kL2(N;R3)= nE

i=1

kyik2

L2(0,`i;R3)!1/2

k(y1,...,ynE)kL2

Hr(N;R3)= nE

i=1

kyik2

Hr(0,`i;R3)!1/2

The norm corresponding to the last term for r= 1 is also used as the norm for H1(N;R6). For

a function y= (y1,...,ynE)∈L2

H1(N;R3) by y0we denote (∂sy1, . . . , ∂synE)∈L2(N;R3).

The results that we prove for the continuous model in Section 3 hold for general geometries,

while the results for the discrete approximation in Section 4 are proved only for stent geometries

with straight struts.

The paper is organized as follows. In Section 2 we present the new formulation of the stent

model. In Section 3 we analyze the new model and give a proof of the inf-sup inequality for

the weak formulation of the continuous infinite dimensional model and in Section 4 we analyze

the discrete model that is obtained after finite element discretization and show a corresponding

inf-sup inequality. Finally in Section 5 we study the dynamical system of the stent movement

under excitation forces. We analyze the properties and present numerical simulation results.

2 New formulation of the model

In this section we reformulate the stent model in a way that enables the proof of a discrete

inf-sup inequality and then, using classical results, appropriate error estimates follow. For this,

we treat all unknowns in the problem explicitly and do not encode them in the function spaces

or the weak formulation. In this way not only the inextensibility and unshearability of the rod

is expressed in the weak formulation, but the continuity of the displacement and infinitesimal

rotation is reflected in the function space of the mixed formulation in H1(N;R6). This leads

to the introduction of new unknowns, displacements and infinitesimal rotations at vertices, and

further, the contact moments and contact forces at the ends of each strut.

Since uand ωare continuous over the whole stent, we introduce as extra variables the

displacements and infinitesimal rotations at the vertices Ui,Ωi,i= 1, . . . , nV, and then form

the vectors

U= [U1,...,UnV]T,Ω= [Ω1,...,ΩnV]T.

Then the kinematic coupling at the vertex jleads to the conditions

ui(`i) = Uj, i ∈J+

j,ui(0) = Uj, i ∈J−

ωi(`i) = Ωj, i ∈J+

j,ωi(0) = Ωj, i ∈J−

j.(2.1)

To express the dynamic coupling conditions, we introduce the contact moments and forces at

the ends of the struts,

+=qi(`i),Qi

−=qi(0),Pi

+=pi(`i),Pi

−=pi(0), i = 1, . . . , nE(2.2)

and define

P±= (P1

±,...,PnE

±),Q±= (Q1

±,...,QnE

±).

Then the dynamic coupling conditions at the vertex jcan be expressed as

i∈J+

+−X

i∈J−

−= 0,X

i∈J+

+−X

i∈J−

−= 0, j = 1, . . . , nE.(2.3)

Equations (1.1)–(1.4), (2.1), (2.2), (2.3), and (1.6) together constitute the stent problem in our

new formulation for which we now derive in detail the weak formulation.

We multiply the ith equation of (1.1) by vi∈H1(0, `i;R3) and that of (1.2) by wi∈

H1(0, `i;R3), add them, integrate over s∈[0, `i], and sum the equations over i. This yields

0 =

i=1 Z`i

∂spi·vi+fi·vi+∂sqi·wi+ti×pi·wids.

After partial integration we obtain

0 =

i=1 Z`i

0−pi·∂svi+fi·vi−qi·∂swi−ti×wi·pids +pi·vi

0+qi·wi

i.e.,

i=1 Z`i

0−pi·(∂svi+ti×wi)−qi·∂swids

+Pi

+·vi(`i)−Pi

−·vi(0) + Qi

+·wi(`i)−Qi

−·wi(0) = −

i=1 Z`i

fi·vids.

(2.4)

In a similar way we multiply (1.3) by ξi∈L2(0, `i;R3) and (1.4) by θi∈L2(0, `i;R3) integrate

over s∈[0, `i] and sum all equations to obtain

0 =

i=1 Z`i

−∂sωi·ξi+Qi(Hi)−1(Qi)Tqi·ξi−(∂sui+ti×ωi)·θids. (2.5)

We also multiply the equations in (2.3) for the jth vertex by Vjand Wjfrom R3, respectively,

and sum the equations over jwhich gives

j=1





X

i∈J+

+−X

i∈J−

−



·Vj+

j=1





X

i∈J+

+−X

i∈J−

−



·Wj= 0.

Since Pi∈J+

jPi

+=Pj

VA+

IP+and Pi∈J−

jPi

−=Pj

VA−

IP−, this equation can be written as

j=1

VA+

IP+−A−

IP−·Vj+

j=1

VA+

IQ+−A−

IQ−·Wj= 0,

and, therefore,

A+

IP+−A−

IP−·V+A+

IQ+−A−

IQ−·W= 0 (2.6)

for all V= [V1,...,VnV]T,W= [W1,...,WnV]T∈R3nV.

Multiplying the equations for the displacements in (2.1) by Θi

+and Θi

−, we obtain

ui(`i)·Θi

+=Uj·Θi

+, i ∈J+

j,ui(0) ·Θi

−=Uj·Θi

−, i ∈J−

Since Pi

E(A+

I)TU=Ujfor i∈J+

j, for Θ±= [Θ1

±,...,ΘnE

±]Twe have

i=1

ui(`i)·Θi

+= (A+

I)TU·Θ+,

i=1

ui(0) ·Θi

−= (A−

I)TU·Θ−,Θ+,Θ−∈R3nE,

and similarly, for the rotations, using the notation Ξ±= [Ξ1

±,...,ΞnE

±]T, we get

i=1

ωi(`i)·Ξi

+= (A+

I)TΩ·Ξ+,

i=1

ωi(0) ·Ξi

−= (A−

I)TΩ·Ξ−,Ξ+,Ξ−∈R3nE.

Thus, we have

i=1

(ui(`i)·Θi

+−ui(0) ·Θi

−)−(A+

I)TU·Θ++ (A−

I)TU·Θ−= 0,Θ+,Θ−∈R3nE,(2.7)

for the displacements and

i=1

(ωi(`i)·Ξi

+−ωi(0) ·Ξi

−)−(A+

I)TΩ·Ξ++ (A−

I)TΩ·Ξ−= 0,Ξ+,Ξ−∈R3nE(2.8)

for the rotations. We multiply the equations (1.6) by αand β, respectively, and summing up,

we obtain

α·ZN

u+β·ZN

ω= 0,α,β∈R3.(2.9)

Subtracting (2.6) from (2.4), we obtain

i=1 Z`i

0−pi·(∂svi+ti×wi)−qi·∂swids

i=1

(Pi

+·vi(`i)−Pi

−·vi(0)) +

i=1

(Qi

+·wi(`i)−Qi

−·wi(0))

−A+

IP+−A−

IP−·V−A+

IQ+−A−

IQ−·W=−

i=1 Z`i

fi·vids,

vi,wi∈H1(0, `i;R3), i = 1, . . . , nE,V,W∈R3nV.

(2.10)

We then add (2.7) and (2.8) to (2.5) and obtain

i=1 Z`i

Qi(Hi)−1(Qi)Tqi·ξi−(∂sui+ti×ωi)·θi−∂sωi·ξids

i=1

(ui(`i)·Θi

+−ui(0) ·Θi

−) +

i=1

(ωi(`i)·Ξi

+−ωi(0) ·Ξi

−)

−((A+

I)TU·Θ+−(A−

I)TU·Θ−)−((A+

I)TΩ·Ξ+−(A−

I)TΩ·Ξ−)=0,

ξi,θi∈L2(0, `i;R3), i = 1, . . . , nE,Θ±,Ξ±∈R3nE.

(2.11)

In [10] and [12] the mixed formulation of the stent model was presented using the space

H1(N;R3) for the displacement vector uand the infinitesimal rotation vector ω. The space

L2(N;R3)×R3×R3for the Lagrange multipliers p,α,β, and the continuity conditions for

displacements and infinitesimal rotations were inherently built into the space H1(N;R3). We

now relax these conditions and consider them as additional equations in the problem and enlarge

the space of unknowns by adding further Lagrange multipliers. The resulting function spaces

are given by

V=L2(N;R3)×L2(N;R3)×R3nE×R3nE×R3nE×R3nE×R3×R3,

M=L2

H1(N;R3)×L2

H1(N;R3)×R3nV×R3nV.

To simplify the notation for the elements of these spaces we introduce

Σ:= (q,p,P+,P−,Q+,Q−,α,β)∈V, φ:= (u,ω,U,Ω)∈M

for the unknowns in the problem and

Γ:= (ξ,θ,Θ+,Θ−,Ξ+,Ξ−,γ,δ)∈V, ψ:= (v,w,V,W)∈M

for the associated test functions.

In this notation, the bilinear forms and the linear functionals that appear in the above

calculations are given by

a:V×V→R, b :V×M→R, f :M→R,

a(Σ,Γ) :=

i=1 Z`i

Qi(Hi)−1(Qi)Tqi·ξids,

b(Σ,ψ) :=

i=1 Z`i

0−pi·(∂svi+ti×wi)−qi·∂swids

i=1

(Pi

+·vi(`i)−Pi

−·vi(0)) +

i=1

(Qi

+·wi(`i)−Qi

−·wi(0))

−A+

IP+−A−

IP−·V−A+

IQ+−A−

IQ−·W+α·ZN

v+β·ZN

f(ψ) := −

i=1 Z`i

fi·vids.

Then the variational formulation (2.10), (2.11) and (2.9) can be expressed as follows.

Determine Σ∈Vand φ∈Msuch that

a(Σ,Γ) + b(Γ,φ) = 0,Γ∈V,

b(Σ,ψ) = f(ψ),ψ∈M. (2.12)

In this way we have obtained that the solution of the stent problem as formulated in (1.1)–(1.6)

satisfies (2.12) and conversely that any solution of (2.12) satisfies (1.1)–(1.6).

In this section we have reformulated the mathematical formulation of the stent model by

including the continuity conditions at the nodes as extra equations and by adding further La-

grange multipliers. In the next sections, we will use this formulation to obtain a discrete inf-sup

inequality and to present a simpler proof of the continuous inf-sup inequality.

3 Properties of the continuous model

In this section we consider the properties of the continuous operator equation (2.12). For the

operator B:V→M0defined by

M0hBΣ,ψiM=b(Σ,ψ),ψ∈M

we have the adjoint operator BT:M→V0(we use the matrix notation to illustrate the

similarity to the discrete case discussed later), which satisfies

VhΣ, BTψiV0=b(Σ,ψ),Σ∈V.

Then Ker BT, the kernel of BT, is defined as a set of vector functions ψ= (v,w,V,W)∈M

such that

b(Σ,ψ)=0,Σ∈V,

so that ψ= (v,w,V,W)∈Ker BTif and only if

∂svi+ti×wi= 0, ∂swi= 0, i = 1, . . . , nE,(3.1)

v=ZN

w= 0,(3.2)

vi(`i) = Pi

E(A+

I)TV,vi(0) = Pi

E(A−

I)TV, i = 1, . . . , nE,(3.3)

wi(`i) = Pi

E(A+

I)TW,wi(0) = Pi

E(A−

I)TW, i = 1, . . . , nE.(3.4)

The conditions (3.3) and (3.4) imply that vand ware continuous on the complete stent, i.e.,

v,w∈H1(N;R3). The conditions in (3.1) imply then that wis constant on the complete stent

and from (3.2) we obtain w= 0. Analogously, from (3.1) we obtain that v= 0, and hence

V=W= 0. Thus we have proved the following lemma.

Lemma 3.1. Ker BT={0}.

As next step we prove that Im BTis closed. For this we derive a kind of Poincar´e inequality

on the graph N, using the notation v(`) = [v1(`1),...,vnE(`nE)]T.

Lemma 3.2. There exists a constant C > 0such that for all v∈L2

H1(N;R3)and all V∈R3nV

the following inequality holds

kvkL2(N;R3)≤Ckv0k2

L2(N;R3)+ZN

v2

+v(0) −(A−

I)TV2

+ v(`)−(A+

I)TV2!1/2

(3.5)

Proof. Suppose the contrary, i.e., for all constants C > 0 there exist vC∈L2

H1(N;R3) and

VC∈R3nVsuch that the opposite inequality holds. Then, for C= 1/k,k∈Nthere exist

sequences vk∈L2

H1(N;R3) and Vk∈R3nVsuch that

kvkkL2(N;R3)= 1,

kv0

kk2

L2(N;R3)+ZN

vk2

+vk(0) −(A−

I)TVk2+vk(`)−(A+

I)TVk2→0,

where, as before, v0

kdenotes the vector of partial derivatives of vk. Thus, taking an appropriate

subsequence (still indexed by k), we have

vk*vweakly in L2(N;R3),(3.6)

k→0 strongly in L2(N;R3),(3.7)

vk→0,(3.8)

vk(0) −(A−

I)TVk→0,(3.9)

vk(`)−(A+

I)TVk→0.(3.10)

It follows that on each strut we have

k*viweakly in L2(0, `i;R3), ∂svi

k→0 strongly in L2(0, `i;R3), i = 1, . . . , nE,

so that viis constant on the ith strut and

k*viweakly in H1(0, `i;R3).

By the Trace Theorem; see e.g. [7, Section 5.5, Theorem 1] we have then

k(0) →vi,vi

k(`i)→vi.

Using (3.9) and (3.10), we have

E(A−

I)TVk→vi,Pi

E(A+

I)TVk→vi, i = 1, . . . , nE.(3.11)

Since in every block row of size 3 the matrices (A−

I)Tand (A+

I)Thave exactly one identity

matrix of size 3 we have that Vkis convergent as well. We denote the limit by V, and have

that its values are given by visuitably organized. Therefore, since AI=A+

I−A−

I, subtracting

the sequences in (3.11) we obtain that

IV= 0.

Since the rank of AT

Iis equal to 3nV−3, the kernel of AT

Iis of dimension 3 by the Rank–Nullity

Theorem, [1]. We easily inspect that Ker AT

Iis spanned by the vectors

(e1,e1,...,e1),(e2,e2,...,e2),(e3,e3,...,e3).

Thus we obtain that all viare equal. Since 0 = RNv=PnE

i=1 `ivi, we obtain that vi= 0 and

hence V= 0, which also implies that v= 0. Thus, since

k(x) = vi

k(0) + Zx

∂svi

k(s)ds,

(vi

k)ktends to 0 strongly in L2(0, `i;R3) for all i= 1, . . . , nE, which is in contradiction to the

unit norm assumption of the sequence, i.e., kvkkL2(N;R3)= 1.

Lemma 3.3. Im BTis closed.

Proof. Consider a convergent sequence in Im BT, i.e., a sequence of the form

∂svi

k+ti×wi

k→pi, ∂swi

k→qi,strongly in L2(0, `i;R3)i= 1, . . . , nE,

vk→α,ZN

wk→β,

k(`i)−Pi

E(A+

I)TVk→Pi

+,vi

k(0) −Pi

E(A−

I)TVk→Pi

−,

k(`i)−Pi

E(A+

I)TWk→Qi

+,wi

k(0) −Pi

E(A−

I)TWk→Qi

−.

(3.12)

Then applying the inequality (3.5) to the sequences wk= (w1

k,...,wnE

k) and Wk= (W1

k,...,WnE

implies that wkis bounded in L2(N,R3). Therefore, wi

kis bounded in H1(0, `i;R3) and hence

there exists a subsequence and a function wi∈H1(0, `i;R3) such that

kl*wiweakly in H1(0, `i;R3), i = 1, . . . , nE.

We collect the limits in w= (w1,...,wnE) and, using again the Trace Theorem, we obtain that

qi=∂swi,−Pi

E(A+

I)TWkl→Qi

+−wi(`i),−Pi

E(A−

I)TWkl→Qi

−−wi(0),

for i= 1, . . . , nEand

β=ZN

Since in each block row of dimension 3 the matrices A±

Ihave exactly one identity matrix, we

obtain that Wklconverges to Wwhich satisfies

−Pi

E(A+

I)TW=Qi

+−wi(`i),−Pi

E(A−

I)TW=Qi

−−wi(0), i = 1, . . . , nE.

By inequality (3.5), there is a unique function wthat satisfies the associated homogeneous

system

∂swi= 0, i = 1, . . . , nE,ZN

w= 0,v(0) −(A−

I)TV=v(`)−(A+

I)TV= 0.

Therefore, wis unique and thus the whole sequences (wk)kand (Wk)kare convergent. Applica-

tion of Lemma 3.2 to wk−wand Wk−Wimplies that wk→wstrongly in L2

H1(N;R3). Once

we have this convergence we apply it to the term ti×wi

kand by the same reasoning identify all

limits related to vkand Vk. We obtain vk→vstrongly in L2

H1(N;R3) and Vk→Vin R3nV

and that

∂svi+ti×wi=pi,vi(`i)−Pi

E(A+

I)TV=Pi

+,vi(0) −Pi

E(A−

I)TV=Pi

−,α=ZN

Thus Σ= (q,p,P+,P−,Q+,Q−,α,β) belongs to Im BT, and hence Im BTis closed.

As a direct consequence of Lemmas 3.1, 3.3, and [4, Proposition 1.2, page 39], we obtain the

following corollary.

Corollary 3.4 (Continuous inf-sup inequality).Consider the variational formulation of the

stent model (2.12). Then there exists a constant kc>0such that

inf

ψ∈Msup

Σ∈V

b(Σ,ψ)

kΣkVkψkM

≥kc.

As our next result we will prove the Ker Bellipticity of the form a, i.e., that there exist

ca>0 such that

a(Σ,Σ)≥cakΣk2

V,Σ∈Ker B.

To obtain this result, we need to restrict the class of networks we consider.

Lemma 3.5. Let the stent geometry be such that

i∈J+

ith edge is straight

αiti−X

i∈J−

ith edge is straight

αiti= 0, j = 1, . . . , nV

implies that αi= 0 for all straight edges i. Then the bilinear form afrom the variational

formulation (2.12) is Ker Belliptic.

Proof. The space Ker Bis given by the set of Σ= (q,p,P+,P−,Q+,Q−,α,β)∈Vsuch that

b(Σ,ψ)=0,ψ∈M.

Thus, from the definition of the form bone has that for all ψ= (v,w,V,W)∈M

i=1 Z`i

0−pi·(∂svi+ti×wi)−qi·∂swids

i=1

(Pi

+·vi(`i)−Pi

−·vi(0)) +

i=1

(Qi

+·wi(`i)−Qi

−·wi(0))

−A+

IP+−A−

IP−·V−A+

IQ+−A−

IQ−·W+α·ZN

v+β·ZN

w= 0.

(3.13)

Since Vand Win RnVare arbitrary, we obtain

A+

IP+−A−

IP−=A+

IQ+−A−

IQ−= 0,(3.14)

which is equivalent to (2.3). These conditions mean that at each vertex the sum of the contact

forces as well as the sum of the contact moments are zero. Then, for γ∈R3, we insert vi=γ,

wi= 0, W=V= 0 as test functions in (3.13) and obtain

i=1

(Pi

+−Pi

−)·γ+α· nE

i=1

`i!γ= 0.

Since from (3.14) nE

i=1

(Pi

+−Pi

−) =

i=1

E(A+

IP+−A−

IP−)=0,

we obtain that α= 0. Now, for fixed iwe insert vi∈C1([0, `i]; R3) with compact support in

h0, `iiin (3.13) and obtain R`i

0pi·∂svi= 0. This implies that piis a constant on each strut.

Inserting a single vi∈C1([0, `i]; R3) in (3.13), we obtain that

−pi·Z`i

∂svids +Pi

+·vi(`i)−Pi

−·vi(0) = 0 (3.15)

and thus pi=Pi

+=Pi

−.

Similarly, for γ∈R3we insert wi=γfor all iin (3.13) and obtain

−

i=1 Z`i

pi×ti·γds +

i=1

(Qi

+−Qi

−)·γ+β· nE

i=1

`i!γ= 0.

As in the case of contact forces we have PnE

i=1(Qi

+−Qi

−) = 0. For the first term we argue as

follows

i=1 Z`i

pi×tids =

i=1

pi×(Φi(`i)−Φi(0)) =

j=1

i∈J+

+−X

i∈J−

−)× Vj= 0,

again by (3.14), where Vjdenotes the jth vertex. Thus we conclude that β= 0. What is left

from (3.13) are the equations for all i∈ {1, . . . , nE}and all wi∈H1(0, `i) given by

−Z`i

pi×ti·wids −Z`i

qi·∂swids +Qi

+·wi(`i)−Qi

−·wi(0) = 0.(3.16)

Defining ˜qi=qi−pi×(Φi(s)−Φi(0)) and inserting this form in (3.16), we obtain

−Z`i

pi×ti·wids −Z`i

˜qi·∂swids −Z`i

pi×(Φi(s)−Φi(0)) ·∂swids

+Qi

+·wi(`i)−Qi

−·wi(0) = 0.

After partial integration in the third term we obtain

−Z`i

˜qi·∂swids + (Qi

+−pi×(Φi(`i)−Φi(0))) ·wi(`i)−Qi

−·wi(0) = 0.(3.17)

As in the case of contact forces, see (3.15), this implies that the ˜qiare constants and that

˜qi=Qi

+−pi×(Φi(`i)−Φi(0)) = Qi

−, i = 1, . . . , nE.

Thus, we have obtained the following characterization of Ker B,

pi=Pi

+=Pi

−,qi=pi×(Φi(s)−Φi(0)) + Qi

−,

+=pi×(Φi(`i)−Φi(0)) + Qi

−,α=β= 0,

with (3.14) being satisfied.

Since Ker Bis of finite dimension and the form ais obviously positive semidefinite, to check

that ais Ker Belliptic, we only have to check that Ker a∩Ker Bis trivial. Assume that

Σ∈Ker a∩Ker B. For elements of Ker a,qi= 0, i = 1, . . . , nEand for elements in Ker Bwe

further have

pi×(Φi(s)−Φi(0)) = Qi

−= 0 s∈[0, `i],Qi

+=α=β= 0, i = 1, . . . , nE.

From the first equation we have that pi=αiti, for some αi∈R,i= 1, . . . , nEif the strut is

straight, otherwise pi= 0. By our assumption on the geometry of the stent this implies that

Ker a∩Ker B={0}. This concludes the proof.

The restriction on the geometry in Lemma 4.4 does not exclude typical examples of stents,

since most of the struts in stents are curved. But even if they were straight, then we can start

from the vertices where only two struts meet and conclude that the associated scalars αifor

these struts should be zero and then continue until we conclude that all coefficients are zero.

The properties of the continuous model that we have shown imply the existence and the

uniqueness of the solution.

Theorem 3.6. There is a unique solution of (2.12).

Proof. Since ais Ker Belliptic by Lemma 3.5 and since the inf −sup inequality holds on V×M

by Lemma 3.1 and Lemma 3.3, the application of [8, Corollary 4.1, page 61] or [4, Theorem 1.1,

page 42] implies the assertion.

Note that, even though we deal with the pure traction problem, because of the introduction

of two additional conditions on total displacement and total infinitesimal rotation in (1.6), we

have a unique solution of (2.12) for all forces. Furthermore, the necessary conditions usual for the

pure traction problem (zero total force and zero total couple) are not necessary any more. The

Lagrange multipliers αand βdeal with that. See [10] for explicit formulas for these multipliers.

4 Properties of the discrete model

In this section we discuss a discrete approximation of the problem (2.12). The derivations in

this section are done for straight struts, i.e., under the assumption that `iti=Φi(`i)−Φi(0) for

all struts. If a strut is curved then it is approximated by a piecewise straight approximation. In

this case we need estimates for two solutions of the stent model for two geometries.

Let us denote by Pm(N) the set of functions on the graph Nwhich are polynomials of degree

m≥0 on each edge of the graph. Note that we do not assume that the functions in Pm(N) are

continuous at vertices. We use the similar notation Pm([0, `]) for the space of polynomials of

degree mon the segment [0, `]. For k, n ∈N0we define the finite dimensional spaces of discrete

approximations

Vk:= Pk(N)3×Pk(N)3×R3nE×R3nE×R3nE×R3nE×R3×R3,

Mn:= Pn(N)3×Pn(N)3×R3nV×R3nV.

We assume that n≥1 and k≥0 and consider the following discrete approximation of (2.12).

Determine Σ∈Vkand φ∈Mnsuch that

a(Σ,Γ) + b(Γ,φ)=0,Γ∈Vk,

b(Σ,ψ) = f(ψ),ψ∈Mn.(4.1)

The form bon Vk×Mndefines the operator Bh:Vk→M0

n, where M0

ndenotes the dual of Mn.

In general, however, Ker Bhis not a subset of Ker B. However, we will show that if n−1≥k

then it is a subset and thus applying Lemma 3.5 gives the following result.

Lemma 4.1. Consider the discrete problem (4.1) and let n−1≥k. Then the bilinear form a

is Ker Bhelliptic.

Proof. As in the continuous case, the elements Σ= (q,p,P+,P−,Q+,Q−,α,β) of Ker Bh

satisfy (3.14) and by the same arguments it follows that α= 0. For fixed iand a test function

vi∈Pn([0, `i]) we obtain the equation

−Z`i

pi·∂svids +Pi

+·vi(`i)−Pi

−·vi(0) = 0.

For constant vi=γ∈R3, we obtain Pi

+=Pi

−and, inserting pi=˜pi+Pi

−implies

−Z`i

˜pi·∂svids = 0,vi∈Pn([0, `i]).

Since n≥k+ 1, the function ˜piis zero and hence pi=Pi

+=Pi

−.

For vi= 0, i= 1, . . . , nEand V=W= 0 in (3.13) we obtain

i=1 Z`i

−pi·ti×wi−qi·∂swids +

i=1

(Qi

+·wi(`i)−Qi

−·wi(0)) + β·ZN

w= 0.(4.2)

Inserting wi=γ∈R3,i= 1, . . . , nE, we obtain

i=1

−γ·Z`i

pi×tids +

i=1

(Qi

+−Qi

−)·γ+β· nE

i=1

`i!γ= 0.

As in the proof of Lemma 3.5, we obtain that β= 0, and thus we are left with the equation

−pi×ti·Z`i

wids −Z`i

qi·∂swids +Qi

+·wi(`i)−Qi

−·wi(0) = 0,wi∈Pn([0, `i]).(4.3)

For k≥1, we insert qi=˜qi+spi×tiinto this equation and obtain

−pi×ti·Z`i

wids −Z`i

˜qi·∂swids −pi×ti·Z`i

s∂swids +Qi

+·wi(`i)−Qi

−·wi(0) = 0.

After partial integration in the third term we obtain

−Z`i

˜qi·∂swids + (Qi

+−`ipi×ti)·wi(`i)−Qi

−·wi(0) = 0.

Since constant functions are contained in Vk, for wi=γwe obtain Qi

+=Qi

−+`ipi×ti. Then

setting ˜qi=˜

˜q+Qi

−, we obtain

Z`i

˜qi·∂swids = 0,wi∈Pn([0, `i]).

As before, since n−1≥k, this implies that ˜

˜q= 0, and hence qi=Qi

−+spi×ti.

For k= 0, qiis constant, so from (4.3) we obtain

−pi×ti·Z`i

wids −qi·(wi(`i)−wi(0)) + Qi

+·wi(`i)−Qi

−·wi(0) = 0.

This implies that

pi×ti= 0,qi=Qi

+=Qi

−

and we have obtained the characterization of Ker Bhgiven by (q,p,P+,P−,Q+,Q−,α,β) that

satisfy (3.14) and

pi=Pi

+=Pi

−,α=β= 0,qi=Qi

−+spi×ti,Qi

+=Qi

−+`ipi×ti.

Additionally, if k= 0, then pi×ti= 0. Thus, Ker Bh= Ker B∩Vk⊆Ker Band hence ais

elliptic on Ker Bhby Lemma 3.5.

Lemma 4.2. Consider the discrete problem (4.1) and let k≥n−1. Then Ker BT

h={0}.

Proof. Ker BT

his defined as a set of ψ= (v,w,V,W)∈Mnsuch that

b(Σ,ψ)=0,Σ= (q,p,P+,P−,Q+,Q−,α,β)∈Vk.

Thus, ψ= (v,w,V,W)∈Ker BT

hif and only if

i=1 Z`i

−pi·(∂svi+ti×wi)−qi·∂swids

i=1

(Pi

+·vi(`i)−Pi

−·vi(0)) +

i=1

(Qi

+·wi(`i)−Qi

−·wi(0))

−A+

IP+−A−

IP−·V−A+

IQ+−A−

IQ−·W+α·ZN

v+β·ZN

w= 0

for all Σ∈Vk. This is equivalent to

i=1 Z`i

0−pi·(∂svi+ti×wi)−qi·∂swids = 0,pi,qi∈Pk([0, `i]),(4.4)

i= 1, . . . , nE,(4.5)

v=ZN

w= 0,(4.6)

vi(`i) = Pi

E(A+

I)TV,vi(0) = Pi

E(A−

I)TV, i = 1, . . . , nE,(4.7)

wi(`i) = Pi

E(A+

I)TW,wi(0) = Pi

E(A−

I)TW, i = 1, . . . , nE.(4.8)

Since k≥n−1, from (4.4) for a test function qiwe obtain that wiis constant for each strut.

The continuity of the infinitesimal rotations at the vertices follows from (4.8). This implies that

wi= const and then (4.6) implies that wi= 0, i.e., wi= 0 for all i= 1, . . . , nEand W= 0.

Analogous arguments using (4.7) imply that vi= 0 and V= 0 as well.

Let n=k+ 1, so that both Lemma 4.1 and Lemma 4.2 apply. Since we are in the finite

dimensional case, clearly Im Bhis closed. Proposition 1.2, page 39 in [4] then implies that

Im Bh= (Ker BT

h)0and then by Lemma 4.2 it follows that Im Bh=M0

n. Furthermore, by

Lemma 4.1, the bilinear form ais Ker Bhelliptic. Then, by the classical theory for finite

dimensional approximations of mixed formulations, e.g. Proposition 2.1 in [4], we obtain the

following existence and uniqueness result for the discretized problem.

Theorem 4.3. Let n=k+ 1. Then problem (4.1) has a unique solution.

Applying the classical results then also we obtain the discrete inf-sup inequality.

Corollary 4.4 (Discrete inf-sup inequality).If n=k+ 1, then there exists a constant kd>0

such that

inf

ψ∈Mn

sup

Σ∈Vk

b(Σ,ψ)

kΣkVkkψkMn

≥kd.

Proof. By Corollary 3.4, the continuous inf-sup inequality holds. By Lemma 3.1 and Lemma 4.2

we have Ker BT

h= Ker BT={0}. Thus Proposition 2.2, page 53 in [4] implies that the as-

sumptions of Proposition 2.8, page 58 in [4] are fulfilled and we obtain the discrete inf-sup

inequality.

Remark 4.1. Note that the constant kdfrom Corollary 4.4 depends on the subspaces Vkand

Mn.

Using Theorem 2.1, page 60 in [4], the discrete inf-sup inequality in Corollary 4.4 and

Lemma 4.1, i.e., the coercivity of the form aon Ker Bh, we obtain error estimates also in the

discrete problem. Introducing analogous notation as in the continuous case,

Σh:= (qh,ph,Ph

+,Ph

−,Qh

+,Qh

−,αh,βh)∈Vk,φh:= (uh,ωh,Uh,Ωh)∈Mn

for the unknowns in the problem and

Γh:= (ξh,θh,Θh

+,Θh

−,Ξh

+,Ξh

−,γh,δh)∈Vk,ψh:= (vh,wh,Vh,Wh)∈Mn

for the test functions we have the following theorem.

Theorem 4.5. Let n=k+ 1 and let (Σ,φ)∈V×Mbe the solution of (2.12) and let

(Σh,φh)∈Vk×Mnbe the solution of (4.1). Then

kΣ−ΣhkV+kφ−φhkM≤cinf

Γh∈Vk

kΣ−ΓhkV+ inf

ψh∈Mn

kφ−ψhkM.

Remark 4.2. The construction of finite elements, as presented, is directly related to the struts

which are described by their prescribed length. To increase the accuracy we can increase the

polynomial degree, assuming that the geometry has been described without error. On the other

hand, we can change the topology of the stent by adding new points on existing struts (and thus

not changing the geometry of the stent) in the original definition of the graph N= (V,E). In

this way we obtain a refined model with transmission conditions of continuity of displacements,

rotations, contact moments and forces at new points. Since at each new vertex only two struts

meet, these transmission conditions are the same coupling conditions (kinematical and dynam-

ical) as for all vertices of the stent. Thus, the resulting weak formulations as in [5] or [10] are

the same for both networks, the original one and the one with added vertices.

The error estimate in Theorem 4.5 can be employed in the finite element method using

interpolation estimates in L2and H1. Note that in contrast to the numerical method in [12],

due to the availability of the discrete inf-sup inequality in the new formulation, we obtain the

error estimate for all variables, including the contact forces. It is a classical result, see e.g. [6],

that for a function ϕ∈Hr(0, `) and its polynomial Lagrange interpolant Πmϕof degree mone

has the estimate

kϕ−ΠmϕkL2(0,`)≤C`min {r,m+1}kϕ(min {r,m+1})kL2(0,`), ϕ ∈Hr(0, `),

kϕ−ΠmϕkH1(0,`)≤C`min {r−1,m}kϕ(min {r,m+1})kL2(0,`), ϕ ∈Hr(0, `).(4.9)

With h:= max{`i, i = 1, . . . , nE}, then combining (4.9) with Theorem 4.5, we obtain the

following error estimate for the finite element method.

Theorem 4.6. Let n=k+ 1,r≥0,r∈N, and let f∈L2

Hr(N;R3). Let (Σ,φ)∈V×Mbe

the solution of (2.12) and (Σh,φh)∈Vk×Mnthe solution of (4.1). Then

kΣ−ΣhkV+kφ−φhkM≤chmin {r+1,k+1}kfkL2

Hr.

Proof. The error of the finite element approximation is estimated by the error of the interpolation

operator. Thus we get

kΣ−ΣhkV+kφ−φhkM≤ckΣ−ΠkΣkV+kφ−ΠnφkM.

Since in this section the struts are assumed to be straight for f∈L2

Hr(N;R3), from the differ-

ential equations we obtain that

p∈L2

Hr+1 (N;R3),q∈L2

Hr+2 (N;R3),ω∈L2

Hr+3 (N;R3),u∈L2

Hr+4 (N;R3).

Thus, we obtain the estimate

kΣ−ΣhkV+kφ−φhkM

≤chmin {r+1,k+1}k(q,p)(min {r+1,k+1})kL2(N;R6)

+hmin {r+2,n}k(u,ω)(min {r+3,n+1})kL2(N;R6)

≤chmin {r+1,k+1}kfkL2

Hmin {r,k}+hmin {r+2,n}kfkL2

Hmin {r,n−2}

≤chmin {r+1,k+1}kfkL2

Hr.

Note that for k= 1 and n= 2 and r≥1 we obtain quadratic convergence in the H1norm

for the displacements and infinitesimal rotations. This is in accordance with the convergence

rate obtained in [12] for the classical formulation.

Having obtained error estimates for the continuous and discrete problem, in the next sub-

section, we consider the properties of the resulting linear system.

4.1 Block structure of the discretization matrix

In the sequel we assume that n=k+ 1. Using the same structure as in the continuous problem,

the discrete problem is given by a linear system Kx=F, where

K=A BT

B0,F=0

F2

with a square matrix Aof size 3(2k+6)nE+6 and a rectangular matrix Bof size (3(2k+4)nE+

6nV)×(3(k+ 6)nE+ 6). Having in mind the evolution problem that we will study in the next

section, we partition these matrices further as

A=A11 0

0 0 ,B=0B32

B41 B42 ,(4.10)

where A11 is a square matrix of size 3(k+1)nE,B32 is a matrix of size 3(k+2)nE×(3(k+5)nE+6),

B41 is of size (3(k+2)nE+6nV)×3(k+1)nEand B42 is of size (3(k+2)nE+6nV)×(3(k+5)nE+6)

associated with the following variables,

dim \unknown q(p,P+,P−,Q+,Q−,α,β)u(ω,U,Ω)F

3(k+ 1)nEA11 0 0 BT

41 0

3(k+ 5)nE+ 6 0 0 BT

32 BT

42 0

3(k+ 2)nE0B32 0 0 F3

3(k+ 4)nE+ 6nVB41 B42 0 0 0

Figure 2: Design of the Palmaz like stent used in simulations

or in more detail,

q p P +P−Q+Q−α β u ω U Ωdimension

ξ⋆0⋆3(k+ 1)nE

θ⋆ ⋆ 3(k+ 1)nE

T+⋆−(A+

I)T3nE

T−⋆(A−

I)T3nE

Ξ+⋆−(A+

I)T3nE

Ξ−⋆(A−

I)T3nE

γ⋆3

δ⋆3

v0⋆ ⋆ ⋆ ⋆ 0 3(k+ 2)nE

ω⋆ ⋆ ⋆ ⋆ ⋆ 3(k+ 2)nE

V−A+

IA−

I3nV

W−A+

IA−

I3nV

4.2 Numerical results

To illustrate our theoretical analysis we test the implementation of the numerical scheme in the

new formulation for a Palmaz type stent as in Figure 2. The radius of the stent is 1.5mm and

the overall length is 1.68cm. There are 144 vertices in the associated graph with 276 straight

edges. All vertices except the boundary ones are junctions of four edges. The cross-sections are

assumed to be square with the side length 0.1mm. The material of the stent is stainless steel

with Young modulus E= 2.1·1011 and Poisson ratio ν= 0.26506. To this structure we apply

the forcing normal to the axis of the of stent, i.e. of the form

f(x) = f(x1)x2e2+x3e3

px2

2+x2

,x= (x1, x2, x3)∈R3,(4.11)

where x1is the axis of the cylinder. As a consequence, the deformation will also posess some

radial symmetry. The problem is a pure traction problem and the applied forces satisfy the

necessary condition. The non-uniqueness of the solution in the problem is fixed using the

Lagrange multipliers αand β.

The solution for the forcing function

f(x1) = 10

105(x1−`/2)2+ 1,(4.12)

is presented in Figure 3; here `is the length of the stent. On the left the solution is projected to

the (x1, x2)–plane, while on the right it is shown from the different perspective. For the forcing

function

f(x1) = 103(x1−`/2)2e3(4.13)

the results are given in Figure 4; again `is the length of the stent.

0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018

−2

−1.5

−1

−0.5

0.5

1.5

x 10−3

00.005 0.01 0.015 0.02

−2

x 10−3

−3

−2

−1

x 10−3

Figure 3: Solution for the radial load from (4.12). On the left the solution is projected to the

(x1, x2)–plane; on the right the solution is shown from a different perspective.

−5 0 5 10 15 20

x 10−3

−3

−2

−1

x 10−3

Figure 4: Solution for the load given in (4.13) projected into (x1, x2)–plane

In the following we present the order of convergence of the finite element method for the

solution of the problem with the quadratic forcing f(x1) = 2.5·107x2

1, the same as in the

numerical scheme presented in [12]. We divide all the edges into 128 smaller rods and solve

the equilibrium problem. The obtained solution we consider as the best possible and use it to

compute the errors, denoted by ”error(i)”, of the approximations for edges split into 2ismaller

struts, i= 1,...,6. We use quadratic finite elements for displacement and infinitesimal rotation,

and linear finite elements for contact forces and couples, i.e. n= 2 and k= 1, see e.g. [6], for

computing the approximations, the L2norm and the H1semi-norm for displacements and the

`1/n norm for unknowns in Rn, i.e., the arithmetic mean of errors, to determine the error

estimates and these to compute the convergence rates via

log error(i+1)

error(i)

log h(i+1)

h(i)

, i = 1,...,5.(4.14)

The obtained convergence rates for u,ω,p,qare in agreement with the analytical estimate

1 1.5 2 2.5 3 3.5 4 4.5 5

1.5

2.5

3.5

4.5

P+

P−

Figure 5: Rate of convergence (4.14) for the L2norm of the contact force pand `1norm of

P+,P−

1 1.5 2 2.5 3 3.5 4 4.5 5

1.5

2.5

3.5

4.5

u in H1

Figure 6: Rate of convergence (4.14) for the H1seminorm and L2norm of the displacement u

and `1norm of U

from Theorem 4.6 for k= 1 and n= 2. In Figure 5 the convergence rates for p,P+and P−are

displayed, while in Figure 6 the convergence rates for the u(L2norm and H1semi-norm) and

Uare plotted. Additionally we present the errors of the remaining unknowns in Table 1.

To compare the numerical scheme for the new formulation with that for the old formulation

in [12] in Table 2, we present the obtained matrix sizes and the computing times (computing

times are determined using difference of ”toc” and ”tic” functions in MATLAB 2010b.) for the

computation on a personal computer with 16GB RAM, 64bit Windows and with INTEL R

Core

i3-7100 [email protected].

The difference between the solution for the displacement u, infinitesimal rotation ωand the

contact force p(the only quantities calculated in the old numerical scheme from [12]) for the

same splitting of edges is at least four digits smaller than the error of the approximation, i.e.

splitting q Q+Q−ωΩ

2 1.8127e-5 3.1945266e-8 3.2871209e-8 8.52118e-4 3.3700235e-5

4 0.4532e-5 0.1995776e-8 0.2024711e-8 1.06516e-4 0.2116646e-5

8 0.1133e-5 0.0124898e-8 0.0125802e-8 0.13314e-4 0.0132275e-5

16 0.0283e-5 0.0007811e-8 0.0007839e-8 0.01664e-4 0.0008260e-5

32 0.0070e-5 0.0000486e-8 0.0000487e-8 0.00208e-4 0.0000514e-5

64 0.0017e-5 0.0000028e-8 0.0000028e-8 0.00025e-4 0.0000030e-5

Table 1: Errors for different splitting of edges (L2errors for qand ω).

new formulation old formulation

splitting # comp. time in ssize of matrix comp. time in ssize of matrix

8 22 105198 2 38958

16 47 211182 42 78702

32 108 423150 152 158190

64 288 847086 629 317166

128 903 1694958 4183 635118

Table 2: Computing times and matrix sizes for the old and new numerical scheme

we obtain the same order of approximation for the same mesh. Further, the size of the matrix

for the new numerics is 2.7 times larger. However, the computing times for the new approach

are smaller when splitting the edges in 32 or more rods.

5 Dynamic modeling of elastic stents

In the previous sections we have considered the static problem for the stent, but for the analysis,

in particular, to study the movement of the stent under the permanent excitation through the

heartbeat, in this section we formulate and analyze the evolution model of the stents.

5.1 Formulation of the space-continuous dynamic model

In order to formulate the dynamic model we start from the evolution equation of curved rods

from [22] and add in the model (1.1)–(1.4) the inertial term to the equilibrium equation (1.1),

which changes to

ρAui

tt =∂spi+fi,

where Ais the area of the cross section and ρis volume density of mass. In the weak formulation

this implies that the term

−

i=1 Z`i

ρAui

tt ·vids

should be added to the left hand side of (2.4) and (2.10). Using the notation of Section 2 and

introducing the bilinear form

c:M×M→R, c(φ,ψ) =

i=1 Z`i

ρAui·vids,

we can formulate the evolution problem of elastic stents as follows:

Determine Σand Γsuch that

a(Σ,Γ) + b(Γ,φ)=0,Γ∈V,

−d2

dt2c(φ,ψ) + b(Σ,ψ) = f(ψ),ψ∈M. (5.1)

5.2 Analysis of the space-discrete model

The dynamical system (5.2) is a finite dimensional linear differential-algebraic equation of second

order, with a nonsingular (but indefinite) stiffness matrix K, and a positive semidefinite but

singular E.

The associated discrete dynamical problem is given by

−E¨

z(t) + Kz(t) = F(t),(5.2)

where the matrix Kand the right hand side Fare as in the static case, Ehas the following

structure partitioned as K

dim \unknown q(p,P+,P−,Q+,Q−,α,β)u(ω,U,Ω)

3(k+ 1)nE0 0 0 0

3(k+ 5)nE+ 6 0 0 0 0

3(k+ 2)nE0 0 M0

3(k+ 4)nE+ 6nV0 0 0 0

and z(t) is the vector of coefficient functions in the finite element basis.

The particular block structure of the DAE allows to analyze the properties of the system,

which are characterized by the spectral properties of the matrix polynomial −λ2E+K.

Lemma 5.1. Consider the DAE (5.2) and the associated pair of matrices (E,K). Then there

exists a nonsingular matrix Vwith the property that

E=VTEV =











,V−1z=











K=VTKV =





0 0 0 0 ˆ

0ˆ

A22 0ˆ

42 0

0 0 ˆ

A33 0 0

0ˆ

B42 0 0 0

B51 0 0 0 0







,VTF=











where ˆ

A33 =ˆ

33,ˆ

B42, and ˆ

B51 are invertible, and ˆ

F4=F3.

Proof. The proof follows by a sequence of congruence transformations, starting from the original

block structure

E=











,K=





A11 0 0 BT

0 0 BT

32 BT

0B32 0 0

B41 B42 0 0







by first compressing

0B32

B41 B42 

via an orthogonal transformation V1from the right to a form

˜

B41 ˜

B42 0

B51 0 0 ,

row-partitioned according to the original row-partitioning, with ˜

B42, and ˜

B51 having full column

rank. Setting ˜

V= diag(V1,I) and applying the transformation with VT

1from the right to the

first two block columns, with V1from the left to the first two block rows, and partitioning

1A11 0

0 0 V1=: 



A11 ˜

A12 ˜

A13

A21 ˜

A22 ˜

A23

A31 ˜

A32 ˜

A33 



accordingly, we obtain a transformed system with

E=˜

VTE˜

V=











,˜

V−1z=











K=˜

VTK˜

V=





A11 ˜

A12 ˜

A13 ˜

41 ˜

A21 ˜

A22 ˜

A23 ˜

42 0

A31 ˜

A32 ˜

A33 0 0

B41 ˜

B42 0 0 0

B51 0 0 0 0







,˜

F:= ˜

VTF=











The property that Ais Ker Belliptic then implies that ˜

A33 is nonsingular. The fact that Kis

invertible and that ˜

B51 has full column rank implies that ˜

B51 is square and nonsingular. Thus,

there exists a nonsingular matrix V2that eliminates the leading 4 blocks in the first block row

and column with ˜

B51 and the leading 2 blocks in the third block column and column with ˜

A33.

Thus, a congruence transformation with V=V2˜

Vyields the asserted structure. The property

that ˜

B42 has full column rank and that ˆ

Kis invertible, then implies that ˆ

B42 =˜

B42 is square

and nonsingular as well.

Remark 5.1. The elimination procedure presented in the proof of Lemma 5.1 is a structured

version of the canonical form construction for differential-algebraic equations, see e.g. [17], which

identifies all explicit or implicit constraints in the system. These lead to restrictions in the initial

values and usually also to further differentiability conditions for the inhomogeneity. Due to the

structure of the equations and the inhomogeneity, the only components of the inhomogeneity

that have to be differentiated in time, are 0, so there are no further requirements on the forcing

function F.

Note further that the same procedure can also be applied in the space-continuous case, by

constructing appropriate projections into subspaces, see e.g. [19]. Since this procedure is rather

technical we do not present this construction here.

Corollary 5.2. Consider the DAE (5.2) transformed as in Lemma 5.1. Then for the general

solution of the transformed system we have ˆ

z1= 0,ˆ

z3= 0,ˆ

z5= 0 ,ˆ

z2is the solution of the

second order DAE ˆ

A22¨

ˆz2=ˆ

42M−1ˆ

B42ˆ

z2+ˆ

42M−1ˆ

F4(5.3)

which exists without any further smoothness requirements for ˆ

F4, and is unique for every con-

sistent initial condition, and finally ˆ

z4=−ˆ

B−T

42 ˆ

A22ˆ

z2.

No initial conditions can be assigned for ˆ

z1,ˆ

z3,ˆ

z5,ˆ

z4, while for ˆ

z2a consistency condition

for ˙

ˆz2in the kernel of ˆ

A22 arises, that depends on the right hand side ˆ

F4.

Proof. This follows directly from the transformed equations. The equations (5.3) form a so

called index-one DAE (see [17]), since ˆ

42M−1ˆ

B42 is positive definite and thus, in particular

invertible in the kernel of ˆ

A22. Projecting into this kernel gives an algebraic equation which has

to hold for the initial condition associated with ˙

ˆz2, while for the remaining components of ˆ

an initial value can can be chosen arbitrarily.

Remark 5.2. The operator pencil associated to the system (5.1) can be studied using the

general theory from [21]. Note that the form e((Σ,φ),(Γ,ψ)) = c(φ,ψ), (Σ,φ),(Γ,ψ)∈

V⊗Mis bounded and symmetric and the form k((Σ,φ),(Γ,ψ)) = a(Σ,Γ) + b(Γ,φ) + b(Σ,ψ),

(Σ,φ),(Γ,ψ)∈V⊗Mis closed, symmetric and semi-bounded from below, see [11]. The form e

defines, in the sense of Kato [16], the bounded, semidefinite and self-adjoint operator E, whereas

the form kdefines the self-adjoint semibounded from below operator K. The operator Kcan

be represented as the formal product K=L∗JL, where Lis a closed operator with a bounded

inverse such that the domains of Land kare equal and Jis the so called fundamental symmetry

(a bounded self adjoint operator such that J2=J). Subsequently it can be shown — using the

special structure of E— that the operator function T(z) = −z2˜

E+˜

K, where ˜

E=L−∗EL−1

and ˜

K=J, is Fredholm operator valued. Furthermore, the operators ˜

Eand ˜

Kare bounded

Hermitian operators and the resolvent set of Tcontains zero by Theorem 3.6. By the results

of [21, Section 1], the pencil Thas finite semi-simple eigenvalues of finite multiplicity. The

construction of an oblique projection onto the reducing subspace associated to the eigenvalue

infinity can be done in a similar way as in Lemma 5.1. The construction is decidedly more

technical and we leave it to a subsequent paper where we will discuss more general second order

systems (e.g. those involving a damping term).

5.3 Numerical results

In this subsection we present some numerical results obtained for the evolution problem. The

time discretization is done using the implicit mid point rule [13] (of convergence order 2) applied

to the first order formulation of the system. At each time step a linear system for the matrix −E+

0.25∆t2Kis solved using backslash in MATLAB. The computations are performed on a personal

computer with 16GB RAM, 64bit Windows and with INTEL R

Core i3-7100 [email protected]. The

presented computing times are determined using the difference of ”toc” and ”tic” functions in

MATLAB 2010b. All simulations are carried out for the following set of parameters:

•elasticity coefficients: µ=E= 1 Pa,

•thickness of the stent struts: 0.0001 m,

•load: radial, as given in (4.11), where

f(x, t) = Fπ

0.003 (x−cvawe(t−t0))(5.4)

with

F(y) = 5·10−8cos(y),if |y|<0.0015

0,else , cvawe = 0.0075, t0= 0.5.

In other words, fis given as traveling wave determined by the function F, where the factor

cvawe denotes the speed of the wave and the term t0asserts the condition f(x, ·) = 0.

•mass density ρ= 2000 kg/m3,

•total time T= 12s.

In Figure 7 the solutions of the problem for the force fin (5.4) at the time-points t∈

{1,2,3,4,5,6}sis plotted.

0 5 10 15 20

x 10−3

−1.5

−1

−0.5

0.5

1.5

x 10−3

0 5 10 15 20

x 10−3

−1.5

−1

−0.5

0.5

1.5

x 10−3

0 5 10 15 20

x 10−3

−1.5

−1

−0.5

0.5

1.5

x 10−3

0 5 10 15 20

x 10−3

−1.5

−1

−0.5

0.5

1.5

x 10−3

0 5 10 15 20

x 10−3

−1.5

−1

−0.5

0.5

1.5

x 10−3

0 5 10 15 20

x 10−3

−1.5

−1

−0.5

0.5

1.5

x 10−3

Figure 7: Solution of the problem projected to the (x1, x2)–plane for loads given in (5.4) and at

times t= 1s, 2s, 3s, 4s, 5s, and 6s.

In the sequel we compare the computing times of two different approaches. In the first

approach we use the MATLAB backslash function to solve the system obtained at every time

step. In the second, we first perform the LDLTdecomposition of the matrix −E+ 0.25∆t2K,

since it is the same in all iterations and then in the time integration we use the obtained LDLT

decomposition to solve the system.

In the Table 5.3 the computing times for two different calculations with different space

discretizations are presented. In the first column the number of splits in the longest edge (strut)

in the stent is displayed. The second column displays the associated number of degrees of

freedom. The remaining columns present the computing times. The results indicate that the

use of the factorization LDLTobviously pays off. Here the time step is equal to 2−4.

In Table 5.3 we compare computing times for different time step sizes with the same final

time T. The total time approximately doubles when lowering ∆twhich is natural, since the

no LDL using LDLT

# splits size of matrix time (s) time for precomputation (s) time for iterations (s)

2012462 476 1.34 110

2125710 554 1.70 203

2252206 793 4.62 392

23105198 1338 23.71 855

Table 3: Computing times for several space discretizations without and with LDLTprecompu-

tation.

number of time steps doubles. However, it is interesting to note that the time for the solution

with LDLTreduces when reducing ∆t. Here every strut is split in 4 smaller struts.

no LDL using LDLT

∆ttime (s) time for precomputation (s) time for iterations (s)

2−3277 3.66 202

2−4554 1.77 392

2−51120 1.52 808

2−62332 1.39 1717

Table 4: Computing times for several values of ∆twithout and with LDLTprecomputation.

We have also computed the errors in the solutions. For the same time step we computed

the errors for different number of strut splits. The errors are presented in Table 5.3. They are

calculated by comparing the solution with the solution for 24strut splits. All errors are presented

in the L2(0, T;L2(N)) norm. In Table 5.3 the errors are given for different ∆tand for fixed

h−1error

200.041110796760794

210.015348501778952

220.003524722458048

230.000224774470258

Table 5: Relative L2(0, T ;L2(N)) errors for different space meshes.

space mesh. The errors are computed with respect to ∆t= 2−7and the same L2(0, T;L2(N))

norm.

∆terror

2−30.038370278764979

2−40.023014149022735

2−50.012051225164378

2−60.003154110021893

Table 6: Relative L2(0, T ;L2(N)) errors for different time steps.

Movies of the dynamic behavior of the stent are presented in the video Appendix.

6 Conclusion

We have presented a new model description for the numerical simulation of elastic stents, which

explicitly displays all constraints. Based on the new formulation an inf-sup inequality for the

finite element discretization is shown and, furthermore, a simplified proof of the inf-sup inequality

for the space continuous problem is presented. Despite an increased number of degrees of

freedom, the new formulation leads to faster simulation time. The presented techniques are

also used to simplify the analysis and numerical solution of the evolution problem describing

the movement of the stent under external forces. Numerical examples illustrate the theoretical

results and show the effectiveness of the new modeling approach.

References

[1] R. B. Bapat,Graphs and matrices, Springer, London; Hindustan Book Agency, New

Delhi, 2014.

[2] C. Beattie, V. Mehrmann, H. Xu, and H. Zwart,Linear port-Hamiltonian descrip-

tor systems, Math. Control Signals Systems, 30 (2018), 17 (27 pages).

[3] D. Boffi, F. Brezzi, and M. Fortin,Mixed Finite Element Methods and Applications,

Springer, Heidelberg, 2013.

[4] F. Brezzi and M. Fortin,Mixed and Hybrid Finite Element Methods, Springer, New

York, 1991.

[5] S. ˇ

Cani´

c and J. Tambaˇ

ca,Cardiovascular stents as PDE nets: 1D vs. 3D, IMA J. Appl.

Math., 77 (2012), pp. 748–770.

[6] P.G. Ciarlet,The Finite Element Method for Elliptic Problems, Society for Industrial

and Applied Mathematics (SIAM), Philadelphia, 2002.

[7] L.C. Evans,Partial Differential Equations, American Mathematical Society, Providence,

1998.

[8] V. Girault and P.-A. Raviart,Finite Element Methods for Navier-Stokes Equations.

Theory and Algorithms, Springer, Berlin, 1986.

[9] G. Griso,Asymptotic behavior of structures made of curved rods, Anal. Appl. (Singap.),

6 (2008), pp. 11–22.

[10] L. Grubiˇ

si´

c, J. Ivekovi´

c, J. Tambaˇ

ca, and B. ˇ

Zugec,Mixed formulation of the one-

dimensional equilibrium model for elastic stents, Rad Hrvat. Akad. Znan. Umjet. Mat.

Znan., 21(532) (2017), pp. 219–240.

[11] L. Grubiˇ

si´

c, V. Kostrykin, K.A. Makarov, and K. Veseli´

c,Representation theo-

rems for indefinite quadratic forms revisited, Mathematika, 59 (2013), pp. 169–189.

[12] L. Grubiˇ

si´

c and J. Tambaˇ

ca,Direct solution method for the equilibrium problem for

elastic stents, Numer. Lin. Algebra Appl., to appear, 2019.

[13] E. Hairer and G. Wanner,Solving Ordinary Differential Equations II: Stiff and

Differential-Algebraic Problems, 2nd revised edition. Springer Verlag, Heidelberg, 1996.

[14] M. Jurak and J. Tambaˇ

ca,Derivation and justification of a curved rod model, Math.

Models Methods Appl. Sci., 9 (1999), pp. 991–1014.

[15] M. Jurak and J. Tambaˇ

ca,Linear curved rod model. General curve, Math. Models

Methods Appl. Sci., 11 (2001), 1237–1252.

[16] T. Kato,Perturbation Theory for Linear Operators, Springer Science & Business Media,

New York, 2013.

[17] P. Kunkel and V. Mehrmann,Differential–Algebraic Equations. Analysis and Numer-

ical Solution, European Mathematical Society (EMS), Z¨urich, 2006.

[18] P. Kunkel, V. Mehrmann, and L. Scholz,Self-adjoint differential-algebraic equa-

tions, Math. Control Signals Systems, 26 (2014), 47–76.

[19] , R. Lamour, R. M¨

arz, and C. Tischendorf,Differential-Algebraic Equations: A

Projector Based Analysis, Springer Science & Business Media, Heidelberg, 2013.

[20] D.S. Mackey, N. Mackey, C. Mehl, and V. Mehrmann,Structured polynomial

eigenvalue problems: good vibrations from good linearizations, SIAM J. Matrix Anal. Appl.,

28 (2006), pp. 1029–1051.

[21] R. Mennicken and M. M¨

oller,Non-Self-Adjoint Boundary Eigenvalue Problems,

North-Holland Publishing Co., Amsterdam, 2003.

[22] J. Tambaˇ

ca,Justification of the dynamic model of curved rods, Asympt. Anal., 31 (2002),

43–68.

[23] J. Tambaˇ

ca, M. Kosor, S. ˇ

Cani´

c, and D. Paniagua,Mathematical modeling of

vascular stents, SIAM J. Appl. Math., 70 (2010), pp. 1922–1952.