Vascular FSI: tissue prestress

$\newcommand{\bm}[1]{\boldsymbol{\mathbf{#1}}} \newcommand{\bvec}[1]{\mathbf{#1}} \newcommand{\bx}{\mathbf{x}} \newcommand{\bd}{\mathbf{d}} \newcommand{\bv}{\dot{\mathbf{d}}} \newcommand{\ba}{\ddot{\mathbf{d}}} \newcommand{\bmw}{\mathbf{w}} \newcommand{\bmv}{\mathbf{v}} \newcommand{\Nh}{N_h} \newcommand{\bU}{\mathbf{U}} \newcommand{\bp}{\bvec{p}} \newcommand{\VV}{\mathbf{V}} \newcommand{\WW}{\mathbf{W}} \newcommand{\QQ}{\mathbf{Q}} \renewcommand{\AA}{\mathbf{A}} \newcommand{\BB}{\mathbf{B}} \newcommand{\CC}{\mathbf{B}} \newcommand{\EE}{\mathbf{E}} \newcommand{\FF}{\mathbf{F}} \newcommand{\FFi}{\mathbf{f}^i} \newcommand{\FFe}{\mathbf{f}^e} \newcommand{\GG}{\mathbf{G}} \newcommand{\CC}{\mathbf{C}} \newcommand{\MM}{\mathbf{M}} \newcommand{\XX}{\mathbf{X}} \newcommand{\KK}{\mathbf{K}} \newcommand{\UU}{\mathbf{U}} \newcommand{\II}{\mathbf{I}} \renewcommand{\SS}{\mathbf{S}} \newcommand{\HH}{\mathbf{H}} \newcommand{\DD}{\mathbf{D}} \newcommand{\NN}{\mathbf{N}} \newcommand{\bmg}{\mathbf{g}} \newcommand{\R}{\mathbb{R}} \newcommand{\JJ}{\mathbf{J}} \newcommand{\RR}{\mathbf{R}} \newcommand{\PP}{\mathbf{P}}$

Blood vessel tissue prestress

In cardiovascular FSI simulations patient-specific geometries of blood vessels are usually reconstructed from in vivo medical imaging such as CT scans. When imaged, the blood vessels are in a state of mechanical stress that puts them in equilibrium with loading coming from the blood flow (i.e. blood pressure and viscous traction). As a result, when modeling blood vessel tissue, a patient specific geometry configuration coming from imaging data may not be used as a stress-free configuration. Taking the in vivo geometry as stress-free is especially not suitable for fluid-structure interaction simulations as this assumption will lead to non-physical large deformations under realistic loading.

To overcome these shortcomings in simulation on patient-specific geometries, techniques that allow for the evaluation of a stress-state that corresponds to a given spatial configuration under known external loads are necessary. Here, we follow the approach proposed by (Hsu and Bazilevs, 2011). To this end, let us recall the weak formulation of the nonlinear elastodynamics equations modeling the solid material in FSI simulations: find $\bd = \bd(t) \in V_D$ such that

where $\Omega^S$ is the solid domain in the reference configuration, $\SS$ is the second Piola–Kirchhoff stress tensor and $\PP = \FF \SS$ is the first Piola–Kirchhoff stress tensor. To account for the prestress state, the variational formulation \eqref{eq:CSM_DynamicWeak} is modified as follows: find $\bd = \bd(t) \in V_D$ such that

$$\begin{equation} \left( \rho \frac{\partial^2 \bd}{\partial t^2} , \bmw \right)_{\Omega^S} + \left( \FF (\SS + \SS_0) , \nabla \bmw \right)_{\Omega^S} = \int_{\Omega} \mathbf{f} \cdot \bmw \, d\Omega + \int_{\Gamma_N} \mathbf{h} \cdot \bmw \, d\Gamma \qquad \forall \bmw \in V, \label{eq:CSM_DynamicWeak2} \end{equation}$$

where $\SS_0$ is an a priori specified prestress tensor. This latter is designed such that in the absence of displacement the blood vessel is in equilibrium with the blood flow forces. This condition leads to the following variational problem, obtained by setting $\bd = \mathbf{0}$ in \eqref{eq:CSM_DynamicWeak2}: find $\SS_0$ such that

where $\Gamma$ is the fluid-structure interface in the reference configuration and $\widetilde{\mathbf{h}}$ is the fluid traction vector. This latter can be obtained, e.g., from a separate rigig-wall blood flow simulation on the reference domain with steady inflow and resistance outflow boundary conditions.

We obtain a particular solution for the state of prestress $\SS_0$ following the iterative procedure detailed in (Hsu and Bazilevs, 2011), which is implemented in the FSI_PrestressSolver function.

A the albegraic level, the prestress term in the structural equations yields the following additional contribute to the solid residual vector:

In particular, the residual vector and the Jacobian matrix become

$$\begin{aligned} \RR^S(\bd_{n+1}^k) &= \frac{1}{\beta \Delta t^2} \MM \bd_{n+1}^k + \FFi (\bd_{n+1}^k) + \RR_P + \JJ_P \bd_{n+1}^k - \FFe_{n+1} - \MM \bm{\xi}_{n},\\[4pt] \JJ^S(\bd_{n+1}^k) &= \frac{1}{\beta \Delta t^2} \MM + \frac{\partial \FFi}{\partial \bd} \left( {\bd_{n+1}^{k}} \right) + \JJ_P, \end{aligned}$$

where

$$(\RR_P)_i = \int_{\Omega^S} \SS_0 : \nabla \bm{\varphi}_i\,d\Omega, \qquad (\JJ_P)_{ij} = \int_{\Omega^S} \nabla \bm{\varphi}_j\SS_0 : \nabla \bm{\varphi}_i\,d\Omega.$$

Note that the vector $\RR_P$ and the matrix $\JJ_P$ do not depend on the unknown displacement. Therefore, they can be assembled once and for all at the beginning of the FSI simulation.

Workflow

The entire procedure for the simulation of an FSI problem with prestress consists of three main steps which are detailed below.

Step 1: CFD simulation with steady inflow

Run a rigid-wall blood flow simulation on the reference domain with steady inflow until convergence to a steady solution

%% Load F mesh
[vertices, boundaries, elements] = msh_to_Mmesh('../mesh/FluidVeryCoarse', dim);

%% Solve Fluid
[U0_Fluid, MESH, DATA] = NSt_Solver(dim, elements, vertices, boundaries, {'P1','P1'}, 'datafile_CFD_Steady', [], 'Figures/PreStressFluid_');

save U0_Fluid U0_Fluid;

In the datafile datafile_CFD_Steady we use the following inflow profile (short ramp to a constant value equal to the flowrate at $t_0 = 0$)

Q_0 = Flow_Rate( 0 );
data.bcDir_t = @(t) Q_0 / 2 * (1 - cos(pi/0.1*t) ) .* (t<0.1) +  Q_0 .* (t>=0.1);

data.bcDir{1,5} = @(x, y, z, t, param)( data.bcDir_t(t)  * N1(1) * 1 + 0.*x.*y); 
data.bcDir{2,5} = @(x, y, z, t, param)( data.bcDir_t(t)  * N1(2) * 1 + 0.*x.*y); 
data.bcDir{3,5} = @(x, y, z, t, param)( data.bcDir_t(t)  * N1(3) * 1 + 0.*x.*y); 

Moreover, we set the ExportFinalLoad field in the datafile in order to compute the traction vector $\widetilde{\mathbf{h}}$ on the FS interface $\Gamma$

data.Output.ExportFinalLoad.computeLoad = true;
data.Output.ExportFinalLoad.flag        = 200;
data.Output.ExportFinalLoad.filename    = 'PrestressFluidLoad.mat';

Step 2: find the prestress tensor $\SS_0$

Launch the FSI_PrestressSolver to find the prestress tensor $\SS_0$ which equilabrates the fluid traction load

[mshS.vertices, mshS.boundaries, mshS.elements, mshS.rings] = msh_to_Mmesh('../mesh/SolidVeryCoarse', dim);
[mshF.vertices, mshF.boundaries, mshF.elements, mshF.rings] = msh_to_Mmesh('../mesh/FluidVeryCoarse', dim);

[R_P, J_P] = FSI_PrestressSolver(dim, mshF, mshS, {'P1','P1'}, 'P1', 'datafile_CFD', 'datafile_CSM_Prestress', [], 'Figures/PreStressSolid_');

Save the residual and Jacobian contributes given by the prestress term in the solid equations, namely $\RR_P$ and $\JJ_P$.

save R_P R_P;
save J_P J_P;

Step 3: run the FSI simulation

Run the FSI simulation with prestresss contribute in the structural equations

%% Solve FSI
FSIt_Solver(dim, mshF, mshS, {'P1','P1'}, 'P1', 'datafile_CFD', 'datafile_CSM', [], 'Figures/FSIsol_');

To this end, the Prestress flag should be set to true in the CSM datafile

data.Prestress = true;

Moreover, initial condition for the fluid velocity and pressure should be loaded from file. To this end, set in the CFD datafile

% initial condition
data.u0{1}    = @(x, y, z, t, param)(load_ICfluid(x, 1));
data.u0{2}    = @(x, y, z, t, param)(load_ICfluid(x, 2));
data.u0{3}    = @(x, y, z, t, param)(load_ICfluid(x, 3));
data.p0       = @(x, y, z, t, param)(load_ICfluid(x, 4));

where

function v = load_ICfluid(x, d)

load U0_Fluid U0_Fluid;
numNodes = length(x);

if d < 4
    v = U0_Fluid(1+(d-1)*numNodes:numNodes*d)';
else
    v = U0_Fluid(1+(d-1)*numNodes:end)';
end
end

References

• M.C. Hsu, and Y. Bazilevs. Blood vessel tissue prestress modeling for vascular fluid–structure interaction simulation. Finite Elements in Analysis and Design, 47(6), 593-599, 2011.
• Y. Bazilevs, K. Takizawa, and T. E. Tezduyar. Computational fluid-structure interaction: methods and applications. John Wiley & Sons, 2013.