Computational Issues for Optimal Shape Design in Hemodynamics

A Fluid-Structure Interaction model is studied for aortic ﬂow, based on Koiter’s shell model for the structure, Navier-Stokes equations for the ﬂuid and transpiration for the coupling. It accounts for wall deformation while yet working on a ﬁxed geometry. The model is established ﬁrst. 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 coeﬃcients 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.

Koiter's Shell Model for Arteries
The following hierarchy of approximations for the displacement d of the aortic wall will be made: • 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 Σ be the shell surface representing the mean position of the blood vessel.
Let n(x) be the normal at x ∈ Σ.Let d(x, t) be the displacement of the wall at x at time t.Normal displacement implies d = ηn.
In [9] it is shown that under such conditions, Koiter's model reduces to the following equation of η on Σ where ρ s is the density and h the thickness of the vessel, T is the pre-stress tensor, 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 −σ s nn .As with all second order wave type equations two conditions must be given at t = 0, for instance leads to the so-called surface pressure model 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, ξ the Poisson coefficient.
Some typical values are (in the metric system MKSA) for a heart beat of one pulsation per second:

Fluid Equations
The Navier-Stokes equations in a moving domain Ω(t) define the velocity u and the pressure p, where ρ f is the density of the fluid and µ its viscosity.Continuity on Σ of fluid and solid velocities implies With the surface pressure model, continuity of normal stresses implies 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 Γ i and Γ o on which we will prescribe pressure and no tangential velocity.If S = Γ i ∪ Γ o , then the boundary Γ is In [10] and many other authors, the matching conditions on Σ are written on the boundary of a fixed reference domain ∂Ω 0 because Koiter's shell model works with a fixed mean surface Σ.
With the notations of [5], assume that the domain of the fluid is Then in Ω t at t = τ , the Navier-Stokes equations are in ALE format 3 Transpiration Conditions for the Fluid

Conservation of Energy
We begin with an important remark on the conservation of energy.
The variational formulation of (3)-divided by ρ s -is, ∀û, p An energy balance is obtained by taking û = u and p = −p, because With transpiration conditions we intend to work on a fixed domain with zero tangential velocity but non zero normal velocity u • n = w.In that case, in order to preserve energy, (6) on a fixed domain Ω needs to be modify into or equivalently into where p = p + 1 2 |u| 2 is the dynamic pressure.Remark 2 Notice that the difference between p and 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 Γ i , Γ o and p Γ 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×n = 0 and shows that we can use several forms for the viscous terms, Hence a variational formulation adapted to the problem is to find u with u×n = 0 and, for all p and all û with û × n = 0

Transpiration
As the wall vessel is {x + ηn : x ∈ Σ} and as, by Taylor, matching the velocities of fluid and structure may be written as On a torus of small radius r and large radius R, at a point of coordinates (R+r cos θ) cos ϕ, (R+r cos θ) cos ϕ, r sin θ), a straightforward calculation shows that So (13) becomes Similarly the normal component of the normal fluid stress tensor is Therefore for a quasi toroidal geometry, for large R, (1) is So, in fine, the domain Ω no longer varies with time but on part of its boundary where a is a non linear function of η.
Remark 3 Notice that η r, i.e. large vessels, allows us to eliminate η and write everything in terms of ∂ t p and u n := u • n.It suffices to differentiate the last equation with respect to t and use the first one and integrate in time,

Variational Formulation and Approximation
Coming back to (12) and using (17): Continuous Problem Find u with u × n = 0 and, for all p and all û with û × n = 0

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, L(u • n) = bu • n.However everything below extends to the full model.
So define Time discrete Problem p(t) = p 0 + bU (t) and we seek where u m+ 1 2 = 1 2 (u m+1 + u m ) and θ = 0 for a semi-explicit lineat but order one scheme or θ = 1 2 for a fully implicit second order scheme in time but nonlinear.

Convergence
A convergence analysis was done in [3]; we recall the results.We denote u δ the linear in time interpolate of Theorem 1 The solution of the time discretized variational problem satisfies Theorem 2 If Ω is simply connected, there is a subsequence (u δ , p δ ) which converges to the continuous problem in L 2 (W) × H −1 (L 2 ) where

Spatial Discretization with Finite Elements
The easiest is to use penalization to enforce u×n = 0 by adding to the boundary integral 1 Σ u m+1 × n • û × 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 Figure 1: An implementation using freefem++ for problem (18) 5 Optimization and Inverse Problems

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 For instance F = |p| 4 will minimize the time averages pressure peak on Σ.

Inverse Problems
Can we recover the structural parameters of the vessel walls from the observation of the pressure?Consider the minimization problem 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.

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 δu, δp.To compute them we linearise the Navier-Stokes equations.These written globally over (0,T) in weak form are, If ûM+1 = 0, r M +1 = 0, it can be rearranged as follows

Adjoint State
To express the variations in terms of δb, we need to introduce an adjoint state v, solution of the following, for all v, q such that v × n = 0 on ∂Ω.
for all v, q such that v × n = 0 on ∂Ω.

Computation of Gradients with Respect to b
Letting v = δu m+1 , q = δp m+1 and summing in m, from 1 to M after multiplication by δt gives, • Display results.

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 ν = 0.01; we chose = ν.The final time is T = 1, the time step is dt = 0.1 and the pressure difference imposed at Γ i (top) and Γ o bottom is 6 cos 2 (π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 Γ i .Consquently the vessel wall becomes fragile as shown by a simulated wall motion by x → x + u m • ndt at every time step, as shown on Figure 5.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).

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.

Figure 2 :
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

Experiment 2 Figure 3 :
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

Figure 4 :
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.

Figure 5 :
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.

Figure 6 :
Figure 6: Left: Optimization criteria Σ×(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.