# Computational Issues for Optimal Shape Design in Hemodynamics

Olivier Pironneau

► To cite this version:

Olivier Pironneau. Computational Issues for Optimal Shape Design in Hemodynamics. Pekka Neittanmaaki. Mathematical modeling and Optimization of complex structures, 40, Springer Verlag, pp.3 - 20, 2015, Computational Methods in Applied Sciences, 10.1007/978-3-319-23564-6_1. hal-01558896

HAL Id: hal-01558896

https://hal.sorbonne-universite.fr/hal-01558896

Submitted on 10 Jul 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Computational Issues for Optimal Shape Design in Hemodynamics

Olivier Pironneau *

March 9, 2015

## Abstract

A Fluid-Structure Interaction model is studied for aortic flow, based on Koiter's shell model for the structure, Navier-Stokes equations for the fluid and transpiration for the coupling. It accounts for wall deformation while yet working on a fixed geometry. The model is established first. Then a numerical approximation is proposed and some tests are given. The model is also used for optimal design of a stent and possible recovery of the arterial wall elastic coefficients by inverse methods.

## Introduction

Hemodynamics, a special branch of computational fluid dynamics, poses many problems of modeling, data acquisition, computation and visualization. However even as of now it is a valuable tool to understand aneurisms, to design stents and heart valves, etc (see for example [13, 6, 12]).

In this paper we shall focus on algorithms for fluid flows with compliant walls like aortic flow, their modelisation, numerical simulation and inverse techniques.

Blood in large vessels like the aorta is Newtonian and flows in a laminar regime with Reynolds number of a few thousands. The Navier-Stokes equation for incompressible fluid is a good model for it.

A blood vessel on the other hand is a complex structure for which linear elasticity is only a first crude approximation and for which the Lamé coefficients do not have a universal value and can vary with individuals.

Nevertheless, like many authors ([11, 9] for instance) we shall use Koiter's linear shell theory.

## 1 Koiter's Shell Model for Arteries

The following hierarchy of approximations for the displacement *d* of the aortic wall will be made:

<footnote>* Sorbonne Universités, UPMC Univ Paris 06, UMR 7598, Laboratoire Jacques-Louis Lions (LJLL), Boite courrier 187, 75252 Paris Cedex 05. France. Olivier.Pironneau@upmc.fr</footnote>
* Small displacement linear elasticity instead of large displacement (needed for the heart).

* No contact inequalities with the surrounding organs.

* Shell model for the mean surface.

* With reference to the mean surface, normal displacement of the walls only.

Let $\Sigma$ be the shell surface representing the mean position of the blood vessel. Let $n(x)$ be the normal at $x \in \Sigma$. Let $d(x, t)$ be the displacement of the wall at $x$ at time $t$. Normal displacement implies $d = \eta n$.

In [9] it is shown that under such conditions, Koiter's model reduces to the following equation of $\eta$ on $\Sigma$

$$ \rho_s h \partial_{tt} \eta - \nabla \cdot (\mathbf{T} \nabla \eta) - \nabla \cdot (\mathbf{C} \nabla \partial_t \eta) + a \partial_t \eta + b \eta = f^s, \quad (1) $$

where $\rho_s$ is the density and $h$ the thickness of the vessel, $\mathbf{T}$ is the pre-stress tensor, $\mathbf{C}$ is a damping term, $a, b$ are viscoelastic terms and $f^s$ the external normal force, i.e. the normal component of the normal stress tensor $-\sigma^s_{nn}$. As with all second order wave type equations two conditions must be given at $t=0$, for instance

$$ \eta|_{t=0} = \eta_0, \quad \partial_t \eta|_{t=0} = \eta'_0. $$

**Remark 1** When $[h, T, C, a] \ll b$, (1) leads to the so-called surface pressure model

$$ -\sigma^s_{nn} = b\eta, \quad \text{with } b = \frac{Eh\pi}{A(1 - \xi^2)}, \quad (2) $$

where $A$ is depends on the geometry of the artery's cross section and equal to the cross section surface when it is circular; $E$ is the Young modulus, $\xi$ the Poisson coefficient.

Some typical values are (in the metric system MKSA) for a heart beat of one pulsation per second:

$$ E = 3\text{MPa}, \quad \xi = 0.3, \quad A = \pi R^2, \quad R = 0.013, \quad h = 0.001, \quad \rho^f = 9.81 \ 10^6, $$

leading to $b = 3.310^7 \text{ms}^{-2}$ and giving displacements in the range of $0.1 \ 10^{-3} \text{m}$ and flow rates around $2 \ 10^{-5} \text{m}^3\text{s}^{-1}$ for aortic flows.

## 2 Fluid Equations

The Navier-Stokes equations in a moving domain $\Omega(t)$ define the velocity $u$ and the pressure $p$,

$$ \rho^f \left( \frac{\partial u}{\partial t} + u \cdot \nabla u \right) + \nabla p - \mu \nabla \cdot (\nabla u + \nabla u^T) = 0, \quad \nabla \cdot u = 0, \quad (3) $$
where $\rho^f$ is the density of the fluid and $\mu$ its viscosity.

Continuity on $\Sigma$ of fluid and solid velocities implies

$$u = \frac{\partial d}{\partial t} := n \frac{\partial \eta}{\partial t}, \text{ on } \Sigma.$$

With the surface pressure model, continuity of normal stresses implies

$$\sigma_{nn}^f := n \cdot (\mu(\nabla u + \nabla u^T) - p)n = -\sigma_{nn}^s := b\eta.$$

Notice that as a consequence of the hypothesis of normal displacements only of the structure, there is no provision to write the continuity of the tangential stresses.

For aortic flow there also an inflow and an outflow boundary $\Gamma_i$ and $\Gamma_o$ on which we will prescribe pressure and no tangential velocity. If $S = \Gamma_i \cup \Gamma_o$, then the boundary $\Gamma$ is

$$\Gamma := \partial\Omega(t) = \Sigma \cup S = \Sigma \cup \Gamma_i \cup \Gamma_o.$$

In [10] and many other authors, the matching conditions on $\Sigma$ are written on the boundary of a fixed reference domain $\partial\Omega_0$ because Koiter's shell model works with a fixed mean surface $\Sigma$.

With the notations of [5], assume that the domain of the fluid is $\Omega_t = A_t(\Omega_0)$ with $A_t: x_0 \to x_t := A_t(x_0)$. Let

$$u_{\tau}(x, t) = u(A_t(A_{\tau}^{-1}(x)), t), \quad \forall x \in \Omega_{\tau}. \qquad (4)$$

Then in $\Omega_t$ at $t = \tau$, the Navier-Stokes equations are in ALE format

$$\rho^f \frac{\partial u_{\tau}}{\partial t} + (u_{\tau} - c_{\tau}) \cdot \nabla u_{\tau} + \nabla p - \mu \nabla \cdot (\nabla u_{\tau} + \nabla u_{\tau}^T) = 0, \\ \nabla \cdot u_{\tau} = 0, \quad \text{with } c_{\tau}(x) = -\frac{\partial A_t(A_{\tau}^{-1}(x))}{\partial t} \Big|_{t=\tau}. \qquad (5)$$

# 3 Transpiration Conditions for the Fluid

## 3.1 Conservation of Energy

We begin with an important remark on the conservation of energy.

The variational formulation of (3)- divided by $\rho^s$ - is, $\forall \hat{u}, \hat{p}$

$$\int_{\Omega(t)} [\hat{u} \cdot (\partial_t u + u \cdot \nabla u) + \nabla p \cdot \hat{u} - \hat{p} \nabla \cdot u + \frac{\nu}{2}(\nabla u + \nabla u^T) : (\nabla \hat{u} + \nabla \hat{u}^T)] = \int_{\Omega(t)} f^s \cdot \hat{u}. \qquad (6)$$

An energy balance is obtained by taking $\hat{u} = u$ and $\hat{p} = -p$,

$$\partial_t \int_{\Omega(t)} \frac{u^2}{2} + \frac{\nu}{2} \int_{\Omega} |\nabla u + \nabla u^T|^2 = \int_{\Omega} f^s \cdot \hat{u} - \int_{\partial\Omega} p u \cdot n, \qquad (7)$$
because

$$ \begin{aligned} \partial_t \int_{\Omega(t)} u \cdot w &= \int_{\Omega(t)} \partial_t (u \cdot w) + \int_{\partial\Omega} v u \cdot w, \\ \int_{\Omega} ((u \nabla u) \cdot u) &= \int_{\partial\Omega} u \cdot n \frac{u^2}{2} = \int_{\partial\Omega} \frac{v}{2} u \cdot u, \end{aligned} \quad (8) $$

when $v = u \cdot n$, the normal speed of $\partial\Omega$.

With transpiration conditions we intend to work on a fixed domain with zero tangential velocity but non zero normal velocity $u \cdot n = w$. In that case, in order to preserve energy, (6) on a fixed domain $\Omega$ needs to be modified into

$$ \begin{aligned} & \int_{\Omega} [\hat{u} \cdot (\partial_t u + u \cdot \nabla u) + \nabla \tilde{p} \cdot \hat{u} - \hat{p} \nabla \cdot u \\ & + \frac{\nu}{2}(\nabla u + \nabla u^T) : (\nabla \hat{u} + \nabla \hat{u}^T)] - \int_{\partial\Omega} \frac{w}{2} u \cdot \hat{u} = \int_{\Omega} f^s \cdot \hat{u}; \end{aligned} \quad (9) $$

or equivalently into

$$ \begin{aligned} & \int_{\Omega} [\hat{u} \cdot (\partial_t u - u \times \nabla \times u) + \nabla \tilde{p} \cdot \hat{u} - \hat{p} \nabla \cdot u \\ & + \frac{\nu}{2}(\nabla u + \nabla u^T) : (\nabla \hat{u} + \nabla \hat{u}^T)] = \int_{\Omega} f^s \cdot \hat{u}, \end{aligned} \quad (10) $$

where $\tilde{p} = p + \frac{1}{2}|u|^2$ is the dynamic pressure.

**Remark 2** Notice that the difference between $p$ and $\tilde{p}$ is second order with respect to the displacement, so exchange one for the other in the shell model is a modification well within the small displacement hypothesis. However it makes a difference on $\Gamma_i$, $\Gamma_o$ and $p_\Gamma$ should be changed accordingly.

From now on we drop the tilde on $p$.

Finally we recall an identity (see [4] for instance) which holds whenever $u \times n = 0$ and shows that we can use several forms for the viscous terms,

$$ \begin{aligned} & \int_{\Omega} [\nabla \times u \cdot \nabla \times v + \nabla \cdot u \nabla \cdot v] = \int_{\Omega} \nabla u : \nabla v \\ & = \int_{\Omega} [\frac{1}{2}(\nabla u + \nabla u^T) : (\nabla v + \nabla v^T) - \nabla \cdot u \nabla \cdot v]. \end{aligned} \quad (11) $$

Hence a variational formulation adapted to the problem is to find $u$ with $u \times n = 0$ and, for all $\hat{p}$ and all $\hat{u}$ with $\hat{u} \times n = 0$

$$ \begin{aligned} & \int_{\Omega} [\hat{u} \cdot (\partial_t u - u \times \nabla \times u) - p \nabla \cdot \hat{u} - \hat{p} \nabla \cdot u + \nu \nabla \times u \cdot \nabla \times v] \\ & + \int_{\partial\Omega} p u \cdot n = \int_{\Omega} f^s \cdot \hat{u}. \end{aligned} \quad (12) $$
## 3.2 Transpiration

As the wall vessel is {$x + \eta n : x \in \Sigma$} and as, by Taylor,

$$u(x + \eta n) = u(x) + \eta \nabla u \cdot n(x) + o(\eta),$$

matching the velocities of fluid and structure may be written as

$$u + \eta \frac{\partial u}{\partial n} = n \frac{\partial \eta}{\partial t} + o(\eta) \text{ on } \Sigma, \quad u \times n = 0. \qquad (13)$$

On a torus of small radius $r$ and large radius $R$, at a point of coordinates $(R+r \cos \theta) \cos \varphi, (R+r \cos \theta) \cos \varphi, r \sin \theta)$, a straightforward calculation shows that

$$u \times n = 0, \quad \nabla \cdot u = 0 \Rightarrow n \cdot \frac{\partial u}{\partial n} = (1 + \frac{r}{R} \cos^2 \theta) \frac{u \cdot n}{r}.$$

So (13) becomes

$$u \cdot n = \partial_t \eta \left(1 + \frac{\eta}{r} \left(1 + \frac{r}{R} \cos^2 \theta\right)\right)^{-1}, \quad u \times n = 0. \qquad (14)$$

Similarly the normal component of the normal fluid stress tensor is

$$\sigma_{nn}^{f} = p + 2\left(1 + \frac{r}{R}\cos^2\theta\right)\frac{\mu}{r}u \cdot n.$$

Therefore for a quasi toroidal geometry, for large $R$, (1) is

$$\rho_s h \partial_{tt} \eta - \nabla \cdot (\mathbf{T} \nabla \eta) - \nabla \cdot (\mathbf{C} \nabla \partial_t \eta) + a \partial_t \eta + b \eta = p + 2 \left(1 + \frac{r}{R} \cos^2 \theta\right) \frac{\mu}{r} \partial_t \eta \left(1 + \frac{\eta}{r} \left(1 + \frac{r}{R} \cos^2 \theta\right)\right)^{-1}. \quad (15)$$

So, in fine, the domain $\Omega$ no longer varies with time but on part of its boundary

$$u \cdot n = \partial_t \eta \left(1 + \frac{\eta}{r} \left(1 + \frac{r}{R} \cos^2 \theta\right)\right)^{-1}, \quad u \times n = 0, \\ \rho_s h \partial_{tt} \eta - \nabla \cdot (\mathbf{T} \nabla \eta) - \nabla \cdot (\mathbf{C} \nabla \partial_t \eta) + a \partial_t \eta + b \eta = p, \qquad (16)$$

where $a$ is a non linear function of $\eta$.

**Remark 3** Notice that $\eta \ll r$, i.e. large vessels, allows us to eliminate $\eta$ and write everything in terms of $\partial_t p$ and $u_n := u \cdot n$. It suffices to differentiate the last equation with respect to $t$ and use the first one and integrate in time,

$$p = p_0 + \mathcal{L}(u \cdot n) := \int_0^t \left( \rho_s h \partial_{tt} u_n - \nabla \cdot (\mathbf{T} \nabla u_n) - \nabla \cdot (\mathbf{C} \nabla \partial_t u_n) + a \partial_t u_n + b u_n \right). \quad (17)$$
# 4 Variational Formulation and Approximation

Coming back to (12) and using (17):

**Continuous Problem** Find $u$ with $u \times n = 0$ and, for all $\hat{p}$ and all $\hat{u}$ with $\hat{u} \times n = 0$

$$ \int_{\Omega} [\hat{u} \cdot (\partial_t u - u \times \nabla \times u) - p \nabla \cdot \hat{u} - \hat{p} \nabla \cdot u + \nu \nabla \times u \cdot \nabla \times v] + \int_{\Sigma} (p_0 + \mathcal{L}(u \cdot n))u \cdot n = - \int_{S} p_{\Gamma} \hat{u} \cdot n $$

## 4.1 Approximation in Time

From now on, for clarity, we consider only the case of the surface pressure model, i.e. $h = T = C = a = 0$, $\mathcal{L}(u \cdot n) = bu \cdot n$. However everything below extends to the full model.

So define

$$ U(t) = \int_{0}^{t} u(s)ds \quad \text{and use the integration rule } U^{m+1} = U^{m} + u^{m+1}dt. $$

and

$$ V = \{u \in H^1(\Omega)^d : u \times n = 0 \text{ on } \partial\Omega\}, \quad Q = L^2(\Omega). $$

**Time discrete Problem** $p(t) = p_0 + bU(t)$ and we seek $u^{m+1} \in V, \hat{p}^{m+1} \in Q$, satisfying for all $\hat{u} \in V, \hat{p} \in Q$,

$$ \int_{\Omega} \left[ \hat{u} \cdot \left( \frac{u^{m+1} - u^m}{\delta t} - u^{m+\frac{1}{2}} \times \nabla \times u^{m+\theta} \right) - p^{m+1} \nabla \cdot \hat{u} - \hat{p} \nabla \cdot u^{m+\frac{1}{2}} + \nu \nabla \times u^{m+\frac{1}{2}} \cdot \nabla \times \hat{u} \right] + \int_{\Sigma} [b\hat{u} \cdot n(u^{m+\frac{1}{2}}\delta t + U^m) \cdot n] = - \int_{S} p_{\Gamma} \hat{u} \cdot n, \quad (18) $$

where $u^{m+\frac{1}{2}} = \frac{1}{2}(u^{m+1} + u^m)$ and $\theta = 0$ for a semi-explicit lineat but order one scheme or $\theta = \frac{1}{2}$ for a fully implicit second order scheme in time but nonlinear.

## 4.2 Convergence

A convergence analysis was done in [3]; we recall the results. We denote $u_\delta$ the linear in time interpolate of $\{u^m\}_1^M$ on $(0, T) = \cup_1^M [(m-1)\delta t, m\delta t]$. For clarity let's assume that $S = \emptyset$.

**Lemma 1** If $\Omega$ is $C^{1,1}$ or polyhedral and $u_0 \in L^2(\Omega)^3$, $p_0 \in H^{1/2}(\Sigma)$, then the weak solution of the continuous problem verifies $u \in L^2(\mathbf{H}^2)$, $\partial_t u \in L^2(\mathbf{L}^2)$, $p \in L^2(\mathbf{H}^1)$, and $u \times n = 0$ in $L^2(L^4(\Sigma))$, $\partial_t p = bu \cdot n$ in $L^2(H^{1/2}(\Sigma))$, $p(0) = p_0$.
**Theorem 1** The solution of the time discretized variational problem satisfies

$$ \|u_{\delta}\|_{L^{\infty}(\mathbf{L}^2)} + \sqrt{\nu} \|u_{\delta}\|_{L^2(\mathbf{H}^1)} + b \|\delta t\|_{L^{\infty}(\mathbf{L}^2(\Sigma))} \\ \qquad \sum_{k=1}^{n+1} u^k \cdot n \|u_{\delta}\|_{L^{\infty}(\mathbf{L}^2(\Sigma))} \\ \qquad \le C \left( \|u_0\|_{0,2,\Omega} + \frac{1}{\sqrt{\nu}} \|p_0\|_{L^2(\Sigma)} \right). $$

**Theorem 2** If $\Omega$ is simply connected, there is a subsequence $(u_{\delta'}, p_{\delta'})$ which converges to the continuous problem in $L^2(\mathbf{W}) \times H^{-1}(L^2)$ where

$$ \mathbf{W} = \{w \in L^2(\Omega) | \nabla \times w \in L^2(\Omega), \nabla \cdot w \in L^2(\Omega), n \times \mathbf{w}_{|_{\Sigma}} = \mathbf{0}\}. $$

## 4.3 Spatial Discretization with Finite Elements

The easiest is to use penalization to enforce $u \times n = 0$ by adding to the boundary integral $\frac{1}{\epsilon} \int_{\Sigma} u^{m+1} \times n \cdot \hat{u} \times n$. Then we may use conforming triangular or tetrahedral elements $P^2$ or $P^1$+bubble for the velocities and $P^1$ for the pressure. A freefem++ implementation (see [8]) is shown on Figure 1.

```c
problem bb([u,v,w,p],[uh,vh,wh,ph], solver=UMFPACK)
    = int3d(th)((u+uh+v+vh+wh)/dt2
                + nu(dx(u)+dx(uh)+dy(u)+dy(uh)+dz(u)+dz(uh) // rot u.rot uh
                + dx(v)+dx(vh)+dy(v)+dy(vh)+dz(v)+dz(vh)
                + dx(w)+dx(wh)+dy(w)+dy(wh)+dz(w)+dz(wh)
                -(dx(u)+dy(v)+dz(w))(dx(uh)+dy(vh)+dz(wh))) // ph div u - p div uh
                -(dx(uh)+dy(vh)+dz(wh))*p - (dx(u)+dy(v)+dz(w))*ph // ph div u - p div uh
                - v*(dx(vold)-dy(vold))*uh - w*(dy(vold)-dz(vold))*vh - u*(dx(vold)-dx(vold))*wh // u x rot u
                + w*(dx(vold)-dx(vold))*uh + u*(dx(vold)-dy(vold))*vh + v*(dy(vold)-dx(vold))*wh // v x rot u
                )
                - int3d(th)(( uold+uh + vold+vh + wold+wh)/dt2 ) // u^m/dt
                + int2d(th,1)(( (u+N.y-v+N.x)*(uh+N.y-vh+N.x) // uxn.uhxn/eps
                + (v+N.z-w+N.y)*(vh+N.z-vh+N.y) // uxn.uhxn/eps
                + (w+N.x-u+N.z)*(wh+N.x-uh+N.z) )/eps)
                + int2d(th,1)(b+dt2*(u+N.x+v+N.y+w+N.z)*(uh+N.x+vh+N.y+wh+N.z)) // b dt u.n uh.n
                + int2d(th,1)(b+(Uold+N.x+Vold+N.y+Wold+N.z)*(uh+N.x+vh+N.y+wh+N.z)) // b n. int_0^t u
                - int2d(th,2)(p0+wh) - int2d(th,3)(p1+vh)
                + on(2,3,u=0) + on(2,v=0) + on(3,w=0);
```

<caption>Figure 1: An implementation using freefem++ for problem (18)</caption>

# 5 Optimization and Inverse Problems

## 5.1 Optimal Stents with the Surface Pressure Model

A stent is a device to reinforce part of a cardiac vessel and/or to change the topology of the flow by its rigidity. This results in a change of the coefficient $b$. So with a first order scheme in time we can consider

$$ \min_{b(x)} J = \int_{\Sigma \times (0,T)} F(p)dxdt \quad : \text{Subject to} $$
$$ \int_{\Omega} [\hat{u} \cdot (\frac{u^{m+1} - u^m}{\delta t} - u^{m+1} \times \nabla \times u^m) - p^{m+1} \nabla \cdot \hat{u} - \hat{p} \nabla \cdot u^{m+1}] \\ + \int_{\Omega} \nu \nabla \times u^{m+1} \cdot \nabla \times \hat{u} + \int_{\Sigma} [b\hat{u} \cdot n(u^{m+\frac{1}{2}}\delta t + U^m) \cdot n] = - \int_{S} p_{\Gamma} \hat{u} \cdot n \\ \forall \hat{u} \in V_h, \hat{p} \in Q_h. \qquad (19) $$

For instance $F = |p|^4$ will minimize the time averages pressure peak on $\Sigma$.

## 5.2 Inverse Problems

Can we recover the structural parameters of the vessel walls from the observation of the pressure?

Consider the minimization problem

$$ \min_{b(x), x \in \Sigma} J(u, p, b) := \frac{1}{2} \int_{\Omega \times (0,T)} (p^m - p_d^m)^2, \qquad (20) $$

subject to (19) or to

$$ \int_{\Omega} [\hat{u} \cdot (\frac{1}{\delta t}(u^{m+1} - u^m(x - u^m(x)\delta t)) - p^{m+1} \nabla \cdot \hat{u} - \hat{p} \nabla \cdot u^{m+1}] \\ + \int_{\Omega} \nu \nabla \times u^{m+1} \cdot \nabla \times \hat{u} + \int_{\Sigma} b(u^{m+1} \delta t + U^m) \cdot \hat{u} = - \int_{\Gamma} p_{\Gamma} \hat{u}_n \\ \forall \hat{u} \in V_h, \hat{p} \in Q_h \text{ with } \hat{u} \times n|_{\Gamma} = 0; \quad U^{m+1} = U^m + u^{m+1} \delta t. \qquad (21) $$

The difference between (19) and (21) is the numerical treatment of the nonlinear term: implicit Euler in the first and Characteristic-Galerkin in the second.

## 5.3 Calculus of Variations

To set up a descent algorithm we must do a sensitivity analysis of the problem. This is done with a “Calculus of Variations”.

When a parameter varies it triggers a variation of $u, p$ which we call $\delta u, \delta p$. To compute them we linearise the Navier-Stokes equations. These written globally over $(0,T)$ in weak form are,

$$ \sum_{0}^{M-1} \delta t \left( \int_{\Omega} [\hat{u}^{m+1} \cdot (\frac{1}{\delta t}(\delta u^{m+1} - \delta u^m(x - u^m(x)\delta t)) - \delta p^{m+1} \nabla \cdot \hat{u}^{m+1} - \hat{p}^{m+1} \nabla \cdot \delta u^{m+1}] \right. \\ + \int_{\Omega} \nu \nabla \times \delta u^{m+1} \cdot \nabla \times \hat{u}^{m+1} + \int_{\Sigma} r^{m+1}(\delta U^{m+1} - \delta U^m - \delta u^{m+1} \delta t) \\ \left. + \int_{\Sigma} (b(\delta u^{m+1} \delta t + \delta U^m) + \delta b(u^{m+1} \delta t + U^m)) \cdot \hat{u}^{m+1} \right) = 0 \\ \forall \hat{u}^m \in V_h, \hat{p}^m \in Q_h, r^m \text{ with } \hat{u}^m \times n|_{\Gamma} = 0. \qquad (22) $$

If $\hat{u}^{M+1} = 0, r^{M+1} = 0$, it can be rearranged as follows

$$ \sum_{0}^{M-1} \delta t \left( \int_{\Omega} \left[ \frac{1}{\delta t} \left( \hat{u}^{m+1} - \hat{u}^{m+2} \left( (x + u^m(x)\delta t) \right) \right) \cdot \delta u^{m+1} - \delta p^{m+1} \nabla \cdot \hat{u}^{m+1} - \hat{p}^{m+1} \nabla \cdot \delta u^{m+1} \right] \right) $$
$$ 
\begin{aligned}
&+ \int_{\Omega} \nu \nabla \times \delta u^{m+1} \cdot \nabla \times \hat{u}^{m+1} + \int_{\Sigma} (r^{m+1} - r^{m+2}) \delta U^{m+1} - r^{m+1} \delta u^{m+1} \delta t \\
&+ \int_{\Sigma} (b(\delta u^{m+1} \delta t + \delta U^m) + \delta b(u^{m+1} \delta t + U^m)) \cdot \hat{u}^{m+1} \bigg) = 0.
\end{aligned} 
\quad (23) $$

## 5.4 Adjoint State

To express the variations in terms of $\delta b$, we need to introduce an adjoint state $v$, solution of the following,

$$ 
\begin{aligned}
& \sum_{0}^{M-1} \delta t \left( \int_{\Omega} \left[ \frac{1}{\delta t} (v^{m+1} - v^{m+2}((x+u^m(x)\delta t)) \cdot \hat{v}^{m+1} - \hat{q}^{m+1} \nabla \cdot v^{m+1} - q^{m+1} \nabla \cdot \hat{u}^{m+1} \right] \right. \\
& \left. + \int_{\Omega} \nu \nabla \times \hat{u}^{m+1} \cdot \nabla \times v^{m+1} + \int_{\Sigma} (r^{m+1} - r^{m+2}) \hat{V}^{m+1} - r^{m+1} \hat{v}^{m+1} \delta t \right) \\
& + \int_{\Sigma} b(\hat{v}^{m+1} \delta t + \hat{V}^m) \cdot v^{m+1} \right) = \sum_{0}^{M-1} \delta t \int_{\Omega} (p^{m+1} - p_d^{m+1}) \hat{q}^{m+1},
\end{aligned} 
\quad (24) $$

for all $\hat{v}, \hat{q}$ such that $\hat{v} \times n = 0$ on $\partial\Omega$. Denote $V^m = r^m \delta t$

$$ 
\begin{aligned}
& \Rightarrow V^{m+1} = V^{m+2} - b v^{m+2} \delta t \\
& \int_{\Omega} \left[ \frac{1}{\delta t} (v^{m+1} - v^{m+2}((x+u^m(x)\delta t)) \cdot \hat{v} - \hat{q} \nabla \cdot v^{m+1} - q^{m+1} \nabla \cdot \hat{v} \right] \\
& + \int_{\Omega} \nu \nabla \times v^{m+1} \cdot \nabla \times \hat{v} + \int_{\Sigma} (b v^{m+1} \delta t - V^{m+1}) \cdot \hat{v} = \int_{\Omega} (p^{m+1} - p_d^{m+1}) \hat{q},
\end{aligned} 
\quad (25) $$

for all $\hat{v}, \hat{q}$ such that $\hat{v} \times n = 0$ on $\partial\Omega$.

## 5.5 Computation of Gradients with Respect to b

Letting $\hat{v} = \delta u^{m+1}$, $\hat{q} = \delta p^{m+1}$ and summing in $m$, from 1 to M after multiplication by $\delta t$ gives,

$$ 
\begin{aligned}
& \sum_{0}^{M-1} \delta t \int_{\Omega} (p^{m+1} - p_d^{m+1}) \delta p^{m+1} \\
&= \sum_{0}^{M-1} \delta t \left( \int_{\Omega} \left[ \frac{1}{\delta t} (v^{m+1} - v^{m+2}((x+u^m(x)\delta t)) \cdot \delta u^{m+1} - \delta p^{m+1} \nabla \cdot v^{m+1} - q^{m+1} \nabla \cdot \delta u^{m+1} \right] + \int_{\Omega} \nu \nabla \times \delta u^{m+1} \cdot \nabla \times v^{m+1} \right. \\
& \left. + \int_{\Sigma} ((b v^{m+1} - r^{m+1}) \delta u^{m+1} \delta t - \delta U^{m+1} (r^{m+2} - r^{m+1} - b v^{m+2})) \right) \\
&= \sum_{0}^{M-1} \delta t \left( \int_{\Omega} \left[ \frac{1}{\delta t} (\delta u^{m+1} - \delta u^m (x - u^m(x)\delta t)) \cdot v^{m+1} \right] \right)
\end{aligned} 
$$
$$ -\delta p^{m+1} \nabla \cdot v^{m+1} - q^{m+1} \nabla \cdot \delta u^{m+1} + \int_{\Omega} \nu \nabla \times \delta u^{m+1} \cdot \nabla \times v^{m+1} \\ + \int_{\Sigma} (bv^{m+1} \delta u^{m+1} \delta t + \delta U^m bv^{m+1} + (\delta U^{m+1} - \delta U^m - \delta u^{m+1} \delta t) r^{m+1}) \\ = - \sum_{0}^{M-1} \delta t \int_{\Sigma} \delta b(u^{m+1} \delta t + U^m) \cdot v^{m+1} = - \int_{\Sigma} \delta b(\delta t \sum_{0}^{M-1} U^m \cdot v^m), \quad (26) $$

because $U^0 = v^M = 0$. To minimize in $H^1$-norm we solve for $g \in H_0^1(\Sigma)$,

$$ \int_{\Sigma} \nabla_s g \cdot \nabla_s w = - \int_{\Sigma} (\delta t \sum_{0}^{M-1} U^m \cdot v^m) w, \forall w \in H_0^1(\Sigma) \\ \Rightarrow \delta J = \int_{\Sigma} \nabla_s g \cdot \nabla_s \delta b \quad (27) $$

## 5.6 Numerical Tests

We take the test case documented in [1]. It is a 2-d problem for the upper part of a symmetric straight vessel. The geometry is the rectangle $(0, L) \times (0, R)$ with $L = 6$ and $R = 0.5$. Pressure is imposed at both end, zero on the right and $p_1 = \frac{1}{2} p_{max}(1 - \cos(2\pi \frac{t}{L}))$ with $p_{max} = 2000$ and $t_{max} = 0.005$.

The mesh is uniform $60 \times 10$. The step size is $\delta t = 210^{-4}$ and there are 60 time steps in this simulation, so $T = 0.012 = 2.4t_{max}$. The $P^2 \times P^1$ element is used for velocity-pressure.

The objective is to see if it is possible to reconstruct $b$ on the upper wall from the pressure in the vessel.

So we first solve the direct problem with $b = b_d := 2.10^5(1 + 6\frac{x}{L}(1 - \frac{x}{L}))$ approximated with the $P^1$ element. We call the computed pressure $\{p_d^m(x)\}_{0}^{M-1}$. Then we solve (20) with 50 iterations of an $H_0^1$-projected gradient method with fixed step size, $\lambda = 10^6$.

### Algorithmic Steps

* Compute $p_d$ by a time loop from 0 to $T$ and store on disk.

* Optimization loop:

1. Compute $u, p$ by a time loop from 0 to $T$ and store on disk $u, p, U$.

2. Compute $v, p$ by a time loop from $T$ down to 0 requiring to read from disk $p_d, u, p, U$.

3. Compute gradient by solving (27).

4. Compute cost function and $\|\partial_x g\|_0^2$.

5. Update $b$ by $b \leftarrow b - \lambda g$.

6. Modify $b$ by $b \leftarrow \max\{\min(b, b_{max}), b_{min}\}$.
<figure><img src="image_2.png" /><figcaption>Figure 2: g (left) and J (right) versus iteration number in log-log scale. Initially $J = 1403$ and after 50 gradient iterations $J = 1.27$ while g decreases from $1.210^{-4}$ to $3.310^{-9}$</figcaption></figure>

• Display results.

We choose $b_{max} = 2.10^5(1 + 12\frac{x}{L}(1 - \frac{x}{L}))$, $b_{min} = 2.10^5(1 + 2\frac{x}{L}(1 - \frac{x}{L}))$.

The results are shown on figures 2, 3, 4.

## 5.7 Preliminary 3D tests

**Experiment 1** This is only a feasibility test with $F = p^4$; The geometry is a quarter of a torus with R=4 and r=1. It is discretized with 1395 vertices and 6120 elements. The number of unknown of the coupled system $[u, p]$ is 23940 with the $P^1$-bubble/$P^1$ element and Crank-Nicolson implicit scheme. The viscosity is $\nu = 0.01$; we chose $\epsilon = \nu$. The final time is $T = 1$, the time step is $dt = 0.1$ and the pressure difference imposed at $\Gamma_i$ (top) and $\Gamma_o$ (bottom) is $6\cos^2(\pi t)$.

The flow is stored on disk at every iteration ready to be reused backward in time for the adjoint equations.

Starting with $b=200$, after 3 iterations of steepest descent with fixed step size, the cost function is decreased from 1200 to 900. But as there is no constraint $b$ is much reduced at the top near $\Gamma_i$. Consequently the vessel wall becomes fragile as shown by a simulated wall motion by $x \to x + \sum u^m \cdot ndt$ at every time step, as shown on Figure 5.

**Experiment 2** The same computations has been made but now $b$ is constrained to be greater than $b_0/2$. A mesh double the size of the previous one has been used, with 191808 degrees of freedom. The initial value of $b$ is $b_0 = 200$. Af-
<figure><img src="image_3.png" /><figcaption>Figure 3: Target $b_d$ (top curve) and computed $b$ after 50 iterations. Initial Pressure map after one iteration (top), final pressure after 50 gradient iterations (middle) and target pressure $p_d$ (bottom). The color scales are linear from -986 to 896 except for $p_0$ which has a range from -680 to 782</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 4: Flow velocity vectors $u$ (middle) and adjoint flow velocity vectors $v$ (bottom) at final time after 50 gradient iterations. The color scales are linear from 0 (safran) to 0.03 (red) for $u$ and 0 (safran) to 2.9 (red) for $v$. The singularity at the top left corner is due to a theoretical incompatibility between the normal velocities at this corner.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 5: Top left: Optimization criteria versus iteration number. Top right: the coefficient b(x) after 3 iterations. Bottom Left: effect of the change of b on the dilatation of the vessel and some iso surfaces of constant pressure. Bottom right: a snap shot of the adjoint pressure and some iso surfaces.</figcaption></figure>
ter 10 iterations, similar to Experiment 1 but with a projected gradient method for the optimization, the results of Figure 6 are found.

<figure><img src="image_6.png" /><figcaption>Figure 6: Left: Optimization criteria $\int_{\Sigma \times (0,T)} p^4$ versus iteration number. Right: the coefficient $b(x)$ after 4 iterations. Right: effect of the change of $b$ on the dilatation of the vessel.</figcaption></figure>

**Experiment 3** Finally we run an identification test of $b$ from the observation of the wall displacement, ideally, $u \cdot n$. However the formulation does not allow it because the extra integral in the adjoint variational formulation is in competition with a similar term from the surface pressure model, so we used $p/b$. For this first test the criteria is

$$J = \int_{\Sigma \times (0,T)} |p - p_d|^2 dxdt,$$

where $p_d$ is obtained from a reference computation (introduction of $b$ in the criteria makes the problem harder) with

$$b = 200 + 100 \cos x \cos y \cos z.$$

The results are shown on figure 7.

Because of the computing cost, we made only an initial study; the target is not reached, but 5 iterations go into the right direction. To do better one would have to used a varying step size gradient method and a better computer (this being done on a macbook pro, takes about 15 min).
<figure><img src="image_7.png" /><figcaption>Figure 7: Left: Optimization criteria $\int_{\Sigma \times (0,T)} (p-p_d)^2$ versus iteration number. Right: the coefficient $b(x)$ after 5 iterations. middle: The target b.</figcaption></figure>

# 6 Conclusion

In this paper we have introduced a reduced fluid structure model based on a transpiration condition and applied it on a problem arising from hemodynamics. We have shown that it has good stability property. In [3] a comparison study is made with full fluid-structure models on moving domains; it is shown to give very similar results.

The greatest advantage of this reduced model is its computational speed and unconditional stability. As inverse problems are important in hemodynamics [2], it could be a good idea to use it. This preliminary study shows that it is indeed feasible.

**Acknowledgement:** *Special thanks to Frédéric Hecht for his help with freefem++ and Marc Thiriet and Suncica Canic for helpful discussions.*

# References

[1] M. BUKACA AND S. CANIC AND R. GLOWINSKI AND J. TAMBACAC AND A. QUAINIA. Fluid-structure interaction in blood flow capturing non-zero longitudinal structure displacement. Journal of Computational Physics 235 (2013) 515-541

[2] C. BERTOGLIO, D. BARBER, N. GADDUM, I. VALVERDE, M. RUTTEN, PH. BEERBAUM, PH. MOIREAU, R. HOSE AND J-F. GERBEAU, Identification of artery wall stiffness: In vitro validation and in vivo results of a data assimilation procedure applied to a 3D fluid-structure interaction model. Journal of Biomechanics 47 (2014) 1027-1034.

[3] T. CHACON REBOLLO, V. GIRault, F. MURAT, O. PIRONNEAU, Analysis of a Simplified Coupled Fluid-Structure Model for Computational Hemodynamics. Submitted to SIAM J Numerical Anal.
[4] M. COSTABEL, M. DAUGE. Singularities of electromagnetic fields in polyhedral domains. Preprint 97-19, Université de Rennes 1 (1997). http://www.maths.univ-rennes1.fr/~dauge/

[5] A. DECOENE AND B. MAURY Moving Meshes with freefem++. J. Numer Math (20)3-4, p195-214(2013).

[6] L. FORMAGGIA, A. QUARTERONI, A. VENEZIANI eds. Cardiovasuclar Mathematics. Springer MS&A series 2009.

[7] V. GIRAULT, R. GLOWINSKI, H. LÓPEZ, J.P. VILA, A boundary multiplier/fictitious domain method for the steady incompressible Navier-Stokes equations, Numer.Math.(2001)88:75-103.

[8] F. HECHT New development in freefem++,nJ. Numer. Math. vol 20, no 3-4, pp 251-265 (2012).

[9] F. NOBILE AND C. VERGANA, an effective fluid-structure interaction formulation for vascular dynamics by generalized robin conditions. SIAM J. Sci. Comp. Vol. 30, No. 2, pp. 731-763 (2008)

[10] F. NOBILE AND C. VERGANA an effective fluid-structure interaction formulation for vascular dynamics by generalized robin conditions. SIAM J. Sci. Comp. Vol.30,No. 2,731-763, 2009

[11] J. TAMBACA, S. CANIC, M. KOSOR, R.D. FISH, D. PANIAGUA. Mechanical Behavior of Fully Expanded Commercially Available Endovascular Coronary Stents. Tex Heart Inst J 2011;38(5):491-501).

[12] M. THIRIET, Biomathematical and Biomechanical Modeling of the Circulatory and Ventilatory Systems. Vol 2: Control of Cell Fate in the Circulatory and Ventilatory Systems. Springer Math& Biological Modeling 2011.

[13] F. USABIAGA, J. BELL, R. BUSCALIONI, A. DONEV, T. FAI, B. GRIFFITH, AND C. PESKIN. Staggered schemes for fluctuating hydrodynamics. Multiscale Model Sim. 10:1369-1408, 2012.