`Unloading`

Module: tutorials.04_EM_tissue.03_unloading.run

Section author: Matthias Gsell <matthias.gsell@medunigraz.at>

Calculate the unloaded reference geometry from a known mesh configuration under loading.

Introduction

Patient-specific models for numerical simulations of the cardiovascular system are mainly based on medical imaging techniques such as X-ray computed tomography (CT) or magnetic resonance imaging (MRI). But at the moment of the medical image acquisition, a physiological pressure load is present. When using the in vivo obtained patient-specific model, this model does not correspond to the unloaded configuration. Thus, we need a suitable algorithm to determine the unloaded configuration from the in vivo model and from the present physiological pressure.

Method

For our purpose, we use a backward displacement method described in [Bols2013] which is based on a fixed-point iteration. Before we describe the method, we make some definitions. We denote by $\Omega(\mathbf{X}, \mathbf{0})$ the stress free reference configuration which is yet unknown. The first argument $\mathbf{X}$ denotes the material coordinates and the second argument $\mathbf{0}$ corresponds to stress ( which is zero ) of the unloaded configuration. A simple forward computation leads to the equilibrium configuration $\Omega(\mathbf{x}, \mathbf{\sigma})$ , with the coordinates of the deformed geometry $\mathbf{x}$ and the second-prder stress tensor $\mathbf{\sigma}$ . This configuration results from a pressure load $p$ applied at the inner surface of the undeformed configuration, i.e.

$p = - \tau \cdot \mathbf{n} = - ( \sigma \, \mathbf{n} ) \cdot \mathbf{n}$

with the outer normal vector $\mathbf{n}$ , and $\tau = 0$ at the outer surface. See figure Fig. 135.

../../_images/unloading_configurations.png

Fig. 135 Unloaded and loaded configuration.

We write

$\Omega(\mathbf{x}, \mathbf{\sigma}) = \mathcal{S}\big(\Omega(\mathbf{X}, \mathbf{0}), p\big)$

where $\mathcal{S}$ is an appropriate forward solver. The deformation is then defined by the mapping $\Phi : \mathbf{X} \mapsto \mathbf{x}$ .

We denote the measured geometry and the measured pressure load by $\mathbf{x}_m$ and $p_m$ respectively. Then the backward problem is as follows.

Find the in vivo configuration $\Omega(\mathbf{x}_m, \mathbf{\sigma}^{\star})$ which is unknown since just $\mathbf{x}_m$ is known, and which is in equilibrium. Therefore, we have to find the corresponding unloaded configuration such that

$\Omega(\mathbf{x}_m, \mathbf{\sigma}^{\star}) = \mathcal{S}\big(\Omega(\mathbf{X}^{\star}, \mathbf{0}), p_m\big)$

where $\mathbf{X}^{\star}$ denotes the unknown unloaded reference geometry which can be written as $\mathbf{X}^{\star} = \Phi^{-1} (\mathbf{x}_m)$ .

Algorithm

The backward displacement method then reads as follows.

$\intertext{\textbf{Input:} Measured state $\Omega_m^r = \Omega(\mathbf{x}_m, \mathbf{0})$ and pressure $p_m$} \\[-1.2cm] & \line(1,0){330} \\[-0.5cm] \intertext{Set initial stress free configuration $\Omega_1^r = \Omega(\mathbf{X}_1, \mathbf{0})$ with $\mathbf{X}_1 = \mathbf{x}_m$} \\[-1cm] & i = 0 \\ & \textbf{while} \; (\text{error} > \varepsilon) \; \textbf{do} \\ & \qquad i = i + 1 \\ & \qquad \Omega_i^t = \Omega(\mathbf{x}_i, \mathbf{\sigma}_i) = \mathcal{S}\big(\Omega_i^r, p_m\big) \\ & \qquad \mathbf{U}_i = \mathbf{x}_i - \mathbf{X}_i \\ & \qquad \mathbf{X}_{i+1} = \mathbf{x}_m - \mathbf{U}_i \\ & \qquad \Omega_{i+1}^r = \Omega(\mathbf{X}_{i+1}, \mathbf{0}) \\ & \qquad \text{error} = error(\Omega_m^r, \Omega_i^t) \\ & \textbf{end do} \\ & \line(1,0){330} \\[-0.5cm] \intertext{\textbf{Output:} Zero pressure geometry $\Omega(\mathbf{X}^{\star}, \mathbf{0})$ with $\mathbf{X}^{\star} = \mathbf{X}_i$} \\[-1.2cm] \intertext{\phantom{\textbf{Output:}} In vivo stress tensor $\mathbf{\sigma}^{\star} = \mathbf{\sigma}_i$ } \\[-1.2cm]$

Usage

To run a simple ring example just run the command

./run.py --meshtype ring --pressure 4.0 --fast --visualize

where the parameter meshtype specifies the mesh and pressure is the pressure to unload from. The visualize flag runs meshalyzer and shows the result. For further information / help on the arguments and flags see below, to get the full help run ./run.py --help.

  -h, --help                             show this help message and exit

  --pressure PRESSURE                    End diastolic pressure to unload from ( in kPa )
                                         (default: 1.5kPa for ring/ellipsoid, 5kPa for cube)

  --np NP                                number of processes

  --meshtype TYPE                        choose the mesh type,
                                         choices = {cube,ring,ellipsoid}

  --dimension DIMENSION                  choose the inner diameter of the mesh cavity, or the
                                         side length in the case of a cube ( in mm )

  --fast                                 set resolution and material (example specific) to
                                         fast-solving values

  --resolution RES                       set target mesh resolution ( in mm )

  --material TYPE                        choose the material model
                                         choices = {mooneyrivlin, linear, neohookean, holzapfelogden,
                                                    demiray, anisotropic, holzapfelarterial,
                                                    stvenantkirchhoff, guccione}

  --parameters PRM=VAL [PRM=VAL ...]     set material model parameters

MATERIAL PARAMETERS:
     mooneyrivlin { kappa, c_1, c_2 }
     linear { E, lmbda, mu, nu }
     neohookean { kappa, c }
     holzapfelogden { kappa, a, a_f, a_fs, a_s, b, b_f, b_fs, b_s }
     demiray { kappa, a, b }
     anisotropic { kappa, a_f, b_f, c }
     holzapfelarterial { kappa, c, k_1, k_2 }
     stvenantkirchhoff { lmbda, mu }
     guccione { kappa, b_f, b_fs, b_t, a }

Examples

In this section we want to present two examples to which the algorithm was applied. The black wireframe shows the initial geometry.

Fig. 136 Three dimensional ring.

Fig. 137 Three dimensional ellipsoid.

References

[Bols2013]

Bols, J. and Degroote, J. and Trachet, B. and Verhegghe, B. and Segers, P. and Vierendeels, J., A computational method to assess the in vivo stresses and unloaded configuration of patient-specific blood vessels (2013), Journal of computational and Applied mathematics, Volume 246