Since the wall shear stress has only L2-regularity in time, we cannot directly use the results fromChapter 3 to prove convergence ofqε to the solution of a suitably averaged equation. While this loss of regularity requires some changes in the convergence proof, we nevertheless use some key results fromChapter 3.
For this, we first prove that the solution v to eq. (4.18) for givenq induces a parametric evolution process as introduced inChapter 3.
Definition 4.3.1. Fort≥s,q∈C(I, Q)andv0∈V define V(t, s;q)v0:=v(t)
where for I := (s, t), v ∈L2(I,D(A))∩H1(I,H),→C(I,V) solveseq. (4.18) with (shifted) initial valuev(s) =v0 and given plaque state q.
Lemma 4.3.2. The set {V(t, s;q);t ≥s, q ∈C(R, Q)} ⊂ C(V,V) is a local-in-time evolution process as defined inChapter 3, i.e. satisfies forq∈C(R, Q) the conditions
• V(t, t;q) =IdV for allt∈R,
• V(t, s;q) =V(t, r;q)V(r, s;q) for allt≥r≥s,
• (t, s, v0)7→V(t, s;q)v0∈V is continuous fort≥s,v0∈V.
Proof. This follows by definition and Lemma 4.2.17using thatL2(I,D(A))∩
H1(I,H),→C(I,V). □
Lemma 4.3.3. There existsα >0such that forq∈C(I, Q)andv1,v2solving eq.(4.18)with given qand initial values v10,v20∈V there holds
k∂t(eαt(v1−v2))kL2(I,H)+keαt(v1−v2)kL2(I,D(A)) ≲kv10−v20kV
and pointwise-in-time exponential stability holds in V fort∈I:
kv1(t)−v2(t)kV≲e−αtkv01−v20kV. Proof. For anyq∈Qthe bilinear form
Vq×Vq 3(v, ϕ)7→(∇v,∇ϕ)q−α(v, ϕ)q
is coercive for0< α < Cq withCq >0 depending only on Poincaré’s constant and hence diam(Ωq). Since diam(Ω)≲diam(Ωq)≲diam(Ω)we can findα >0 independent ofqsuch that the bilinear form is coercive for allq∈Q.
Analogous to the arguments for the bilinear form (4.15) the form above induces a maximal regular operator A˜q,α = ˜Aq−αIdHq on Hq with domain D( ˜Aq)and operator norm equivalent tok · kD( ˜Aq)for allq∈Q, with analogous properties forAq,α:=Pq◦A˜q,α◦Pq−1=Aq−αIdH. Forq∈C(I, Q),t7→Aq(t),α is maximal regular, so that the solution to
∂tw(t) +Aq(t)w(t)−αw(t) = 0, w(0) =w0∈V (4.22) obeys the bound
k∂twkL2(I,H)+kwkL2(I,D(A))≲kw0kV.
Noww:=eαt(v1(t)−v2(t))satisfieseq. (4.22), which follows by subtraction of the equations forv1andv2multiplied with eαt. Thus the maximal regularity estimate forwimplies the first claimed bounds. The embeddingL2(I,D(A))∩ H1(I,H),→C(I,V)results in the pointwise-in-time estimates. □ Corollary 4.3.4. The assumptions A1–A3 fromChapter 3are satisfied, i.e.
A1. There is α >0 such that
kV(t, s;q)v0−V(t, s;q)˜v0kV ≲e−α(t−s)kv0−v˜0kV. A2. There is a monotone function λ:R≥0→R≥0 such that
kV(t, s;q)v0−V(t, s; ˜q)v0kV ≲kq−q˜kC([s,t],Q) kv0kV+λ(t−s) fors < t,v0∈V andq,q˜∈C([s, t], Q).
A3. For eachq∈Qthe process is1-periodic, i.e.
V(t+ 1, s+ 1;q) =V(t, s;q) for anys < t.
Proof. Exponential stability, pointwise inV, was verified inLemma 4.3.3, the Lipschitz continuity inLemma 4.2.18and periodicity follows since the datapio
is1-periodic. □
Lemma 4.3.5. For each q∈ Q there is a unique vπ(·;q) ∈L2(0,1;D(A))∩ H1(0,1;H)solving the time-periodic equation
∂svπ(s;q) +Aqvπ(s;q) =f(s;q), vπ(0, q) =vπ(1, q).
The map
Q3q7→vπ(·;q)∈L2(0,1;D(A))∩H1(0,1;H) is Lipschitz continuous.
Proof. By Lemma 3.3.2a unique periodic solution exists. The Lipschitz con-tinuity in L2(0,1;D(A))∩H1(0,1;H) follows fromLemmas 4.3.3 and 4.2.18 together with the Lipschitz continuity ofQ3q7→vπ(0;q)∈V also proven in
Lemma 3.3.2. □
Let us remind of our notationτ∈I := (0,T)forT >0 to denote the slow timescale. For some fixed0< ε1it is related to the fastε-timescale, where we write t ∈I := (0, T)for T > 0, by τ =εt, T =εT. The limit system is naturally formulated in the slow timescale, but we return to the fast scale for the convergence analysis.
Problem 4.3.6 (Limit system). Given T >0, q0 ∈Q, find q0 ∈C0,1(I, Q) with I:= (0,T)such that forτ∈I there holds
q′0(τ) =g(q0(τ)), q0(0) =q0 where the averaged right-hand side g:Q→Rn is defined by
g(q) :=
Z 1 0
g(q, vπ(s;q))ds (4.23) andvπ(·;q) forq∈Qdenotes the time periodic solution from Lemma 4.3.5.
Lemma 4.3.7. The functiong defined in eq.(4.23) is Lipschitz continuous.
Proof. Forq1, q2∈Qthere holds for i∈ {1,2}:
|g(q1)−g(q2)|
≲ Z 1
0
|q1−q2|(1 +kvπ(s;qi)kD(A)) +kvπ(s;q1)−vπ(s;q2)kD(A)ds
≲|q1−q2|(1 +kvπ(qi)kL2(0,1;D(A))) +kvπ(q1)−vπ(q2)kL2(0,1;D(A))
≲|q1−q2|
using the Lipschitz-continuity ofg,ineq. (4.9), and the Lipschitz continuity of q7→vπ(q)∈L2(0,1;D(A))fromLemma 4.3.5. □
Lemma 4.3.8. For every q0∈Qthere exists T >0 such thatProblem 4.3.6 has a unique solution q0 ∈ C1,1(I, Q) where I := (0,T). For the maximum time of existence eitherT =∞orq0(T)∈∂Q.
Proof. This follows from standard theory of ordinary differential equations due to the Lipschitz continuity and boundedness ofg fromLemma 4.3.7. □
Theorem 4.3.9. Forq0∈Qandv0∈V there existsT >0such that for any 0< ε1 the solution (qε, vε)toProblem 4.2.16 exists onI:= (0, ε−1T), the solution q0 toProblem 4.3.6exists on I= (0,T)and there holds
kqε−q0kC(I,Q)≲ε.
Proof. By Corollary 4.2.20 there exists T >0 such that for any ε >0 Prob-lem 4.2.16 has a solution (qε, vε) on I := (0, T) with T := Tε := ε−1T. By Lemma 4.3.8the same holds true forq0withT smaller if necessary. Fort∈I there holds
|qε(t)−q0(t)|2
≤ε2| Z t
0
g(qε, vε)−g(q0)ds|2
≤ε2| Z t
0
g(qε, vε)−g(q0, vε)ds|2+ε2| Z t
0
g(q0, vε)−g(q0)ds|2. (4.24)
For the first term we use the Lipschitz continuity ofg,ineq. (4.9), to estimate
which together withT ≲ε−1 implies forineq. (4.24)that
|qε(t)−q0(t)|2≲ε Our goal is to apply Gronwall’s inequality to ineq. (4.25), which yields the claimed estimate if all terms on the right are either of the form of the first term or of order ε2. We abbreviatev0π(q) :=vπ(0;q)and ˜vε(t) :=v(0, t;qε)vπ0(q0). and the exponential stability estimate fromLemma 4.3.3 implies
Z t
For the second term on the right ofineq. (4.26)we proceed as inChapter 3and split[0, t]into btcintervals of period length1 and a remainder:
|
The last term is bounded independent ofεsinceg, and henceg, is. We abbre-viate e.g. qεi−1:=qε(i−1) and split the remaining integrals over full periods
as follows: The integral inline (4.27e)is zero by definition ofgsince we integrate over one period. Since we sum overbtc≲ε−1 periods we requirelines (4.27b),(4.27d) and (4.27f) to be of order ε, line (4.27c) will be estimated differently. For line (4.27b)we use Lipschitz continuity of g, such that it remains to estimate
Z i argument in Lemma 3.4.4, noting that not the specific form of qε but only a Lipschitz constant of orderεis necessary there. For the second term we have
Z i theO(ε)-estimate follows sinceg is Lipschitz continuous byLemma 4.3.7. For
the last term online (4.27c)we have by similar techniques as before Z i
i−1
|g(q0, vπ(qεi−1))−g(q0, vπ(q0i−1))|ds≲ Z i
i−1
kvπ(qiε−1)−vπ(q0i−1)kD(A)ds
≲kvπ(qεi−1)−vπ(qi0−1)kL2(0,1;D(A))≲|qiε−1−q0i−1|. Using the Lipschitz continuity ofqεandq0 we have
|qεi−1−qi0−1| ≤ Z i
i−1
|qiε−1−qε|+|qε−q0|+|q0i−1−q0|ds≲ε+ Z i
i−1
|qε−q0|ds.
Applying these estimates forlines (4.27b)to(4.27f)inineq. (4.25), we conclude
|qε(t)−q0(t)|2
≲ε Z t
0
|qε(s)−q0(s)|2ds+ε2+ε2
⌊
t⌋
X
i=1
Z i i−1
|qε(s)−q0(s)|ds
2
≲ε Z t
0
|qε(s)−q0(s)|2ds+ε2
and thus by Gronwall’s inequality the claimed estimate follows:
|qε(t)−q0(t)|2≲ε2eCεt ≲ε2eCT. □
Numerics for the Simplified Plaque Model
This chapter is concerned with the numerical solution of the problem described in Chapter 4, i.e. a fluid equation with time-periodic boundary conditions on a slowly moving domain with movement governed by an ordinary differential equation depending on the wall shear stress. In contrast to the theory we will mainly study the Navier–Stokes equation, assuming that analogous results for the multiscale convergence hold. We emphasize that the interest lies in the evolution of the growth, rather than the blood flow for which a simple reconstruction is available.
The motivation for the theory was that the limit equation is much more efficient to solve numerically than the original fast-slow system, at least for ε
“sufficiently small”. The numerical solution of the limit system is non-standard since each evaluation of the right-hand side of the limit equation forq0requires the solution of the time-periodic solution to the Navier–Stokes equation on the domain corresponding to q0. This chapter consists of a description of the problem and the discretizations used for both fast-slow and limit systems, an error analysis for the temporally semi-discrete limit equation, the numerical convergence results and an improvement, both theoretically and numerically, of an existing method to solve the time-periodic fluid equation.
5.1 Problem Description
For computational simplicity we only study the two-dimensional geometry from Example 4.1.2fromChapter 4, see alsoFigure 5.1. To summarize, this domain is a single-channel geometry with dimensionless diameter 1 and length L :=
5. At the bottom-center a stenosis may develop with base width L0 := 3.5, prescribed bump-shape and height given by the slow variable q. The fluid flows from the inflow boundary ∂iΩ on the left to the outflow boundary∂oΩ on the right, the remaining boundary is the fixed wall∂wΩq. Only the region Γq⊂∂wΩq of the bump is considered damaged and hence permeable, which in our model means that the growth rate depends on the wall shear stress only in this region.
Some physiological parameters such as Reynolds and Strouhal number were taken from the non-dimensionalization of the model by Yang et. al. [Yan+15]
93
q
∂wΩq
3.5 1
∂iΩq ∂oΩq
Γq⊂∂wΩq
5
Figure 5.1: Sketch of the problem domainΩq for someq≈0.5with notation for boundaries and geometric dimensions. The reference domain Ωˆ is a rectangle with the same dimensions. This figure is a repe-tition ofFigure 4.1.
from Section 2.3, but the situation treated here must be viewed as a purely theoretical test case to demonstrate the numerical implementation and accel-eration of the proposed multiscale scheme compared to the direct numerical simulation of the original fast-slow problem. Apart from the geometrical sim-plifications we investigate an exemplary flow profile and wall shear stress to growth rate relation.
We formulate all equations on the unknown moving domain. The prob-lem is then to find fluid velocity vε, pressure pε and growth qε as solution to the following equations. The fluid velocity and pressure is solution to the (dimensionless) Navier–Stokes equation on a moving domain given by
St∂tvε+ (vε· ∇)vε−Re−1∆vε+∇pε= 0, divvε= 0,
)
inΩIqε (5.1a) with Strouhal numberSt= 0.016and Reynolds numberRe= 500in the slowly evolving space-time domain
ΩIqε :={(t, x)|t∈I, x∈Ωqε(t)} (5.1b) with plaque state qε given by
qε′ =εg(qε, vε) (5.1c)
where
g(q, v) :=γ(1−q) Z
Γq
1 + |σWS(v)| σWS0
−1
do
following the example from eq. (4.3) where σWS(v) is the wall shear stress corresponding to the velocity fieldv. The parameterσWS0 := 0.004was chosen such that the integrand over Γq has non-trivial behavior, whereas γ := 0.1 rescales the slow time such that interesting behavior happens for τ ∈ [0,10]
and effectively acts as a rescaling ofε.
The initial values chosen for our computations areqε(0) := 0andvε(0) := 0 in Ω0. From a modelling perspective it would be more appropriate to set vε(0) :=vπ(0; 0)under the assumption that the blood flow prior to our model’s initial time is undisturbed and hence periodic. Such an initial value is inves-tigated for a single test case below, but it requires knowledge of the initial
periodic flow, which is unknown in general1. The method we present to solve the time-periodic Navier–Stokes problem in our limit equation could be used to determine this initial value, but this would convolute the presentation and com-putations. This nevertheless highlights that a time-periodic problem should be solved for realistic models even if the fast-slow problem is studied. We remind that the limit system is independent of the choice of initial fluid velocity, so the simplificationvε(0) := 0mainly affects the boundary layer behavior for the fast variablevε.
On in- and outflow boundaries, a1-periodic pressure drop is prescribed, Re−1∂nvε−pεn=−Pion onI×∂ioΩ (5.1d) with temporally 1-period pressure drop functionPio described below. For the wall boundary homogeneous Dirichlet values are used for simplicity
vε(t, x) = 0 fort∈I andx∈∂wΩqε(t). (5.1e) The pressure drop function is graphed in Figure 5.2and given by
Pio(t, x) = (p0
2 1−cos(2π(2t−t2))
forx∈∂iΩ,
0 forx∈∂oΩ (5.2)
for t ∈ [0,1], periodically extended to R. The parameter p0 := 0.16 was chosen such that the induced periodic velocity field has a maximal speed of approximately 1 in the initial geometry while the solution is still stable, see the remarks about the pressure boundary conditions below. The shape of the pressure drop function is intended to only roughly approximate systole (heart contraction) and diastole (refilling of the heart), including the asymmetry of this process.
Translating the multiscale result fromChapter 4 to the present situation for the Navier–Stokes equation we assume that the limit asε→0is given by the solution to
q′0(τ) =g0(q0(τ)) (5.3a) with initial valueq0(0) := 0and right-hand side
g0(ˆq) :=γ(1−q)ˆ time-periodic Navier–Stokes equation on a cylindrical domain
St∂svπ(ˆq) + (vπ(ˆq)· ∇)vπ(ˆq)−Re−1∆vπ(ˆq) +∇pπ(ˆq) = 0, divvπ(ˆq) = 0,
)
in (0,1)×Ωqˆ, (5.3c) omitting the temporal variable s ∈ (0,1) here, i.e. writing vπ(ˆq) instead of vπ(s; ˆq). This solution is temporally 1-periodic,
vπ(0; ˆq) =vπ(1; ˆq) in Ωqˆ, (5.3d)
1For the straight tube geometry one could use (a spatial projection of) the Womersley flow as an exact solution [Wom55], but such a solution is also only given as a Fourier series depending on the Fourier decomposition ofpio.
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
t
inflow pressure
Figure 5.2: Plot of the pressure (drop) at the inflow boundary∂iΩover one period, seeeq. (5.2). A zero pressure is prescribed on the outflow boundary∂oΩ.
and has the same boundary conditions as the fast-slow problem, i.e.
Re−1∂nvπ(ˆq)−pπ(ˆq)n=−Pion on(0,1)×∂ioΩ, (5.3e) vπ(ˆq) = 0 on(0,1)×∂wΩqˆ. (5.3f)
We refer to [GK16] for an overview of the theory for the time-periodic Navier–Stokes equation, in particular the existence and uniqueness of weak time-periodic solutions with homogeneous Dirichlet boundary conditions for sufficiently regular right-hand side and domain, assuming the smallness of the data ford= 3. Corresponding results for the boundary conditions and domain studied here are not known to the author. In fact, the use of the do-nothing boundary condition on an outflow boundary, as ineq. (5.3e), is known to cause stability issues due to the possibility of uncontrollable backflow, i.e. inflow at nominal outflow boundaries. As a remedy, directional do-nothing conditions have been proposed which add a stabilization only affecting such backflow [BM14] among other approaches [Ber+18]. On inflow boundaries such modifi-cations drastically change the flow profile and are hence unsuitable. Instead, tame problem data and geometry were selected such that no stabilization is necessary.
The behavior of the solutions qε and q0 are plotted in Figure 5.3. These solutions were obtained numerically using the discretization discussed below.
All simulations were stopped after the growth variable exceeded a80%-diameter stenosis, i.e. q > qmax := 0.8. In the algorithmic descriptions of the time-stepping schemes this condition is replaced with a fixed number of time steps for simplicity. All errors between solutions in the following were computed on the largest shared time interval.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
0 1 2 3 4 5 6 7 8 9 10
slow time τ q0
q1 q0.1 q0.01
0.774 0.775 0.776 0.777 0.778
8.56 8.57 8.58 8.59 8.6
Figure 5.3: Behavior of the solutions qε for different values of ε and the limit solution q0. The simulation is stopped after q exceeded qmax:= 0.8(dashed line), the limit system crosses this threshold significantly due to the large step size employed. For smallε, qε
and q0 are indistinguishable on this scale, hence a zoomed view is provided.