The Variational Multiscale Formulation for the Fully-Implicit Log-Morphology Equation as a Tensor-Based Blood Damage Model

TL;DR

This paper introduces a variational multiscale finite element formulation for a tensor-based blood damage model that enhances numerical stability and accuracy, especially in complex blood flow simulations involving ventricular assist devices.

Contribution

The paper develops a novel VMS stabilization method for the log-morph equation, improving numerical behavior over existing GLS and SUPG methods in blood damage modeling.

Findings

VMS stabilization outperforms GLS in 2D stirrer tests.

VMS method shows clear advantages over SUPG in VAD simulations.

Logarithmic tensor description prevents unphysical eigenvalues.

Abstract

We derive a variational multiscale (VMS) finite element formulation for a viscoelastic, tensor-based blood damage model. The tensor equation is numerically stabilized by a logarithmic shape tensor description that prevents unphysical, negative eigenvalues. The resulting VMS stabilization terms for this so-called log-morph equation are presented together with their special numerical treatment. Results for a 2D rotating stirrer test case obtained from log-morph simulations with both SUPG and VMS stabilization show significantly improved numerical behavior if compared with Galerkin/least squares (GLS) stabilized untransformed morphology simulation results. The newly proposed method is also successfully applied to a state-of-the-art centrifugal ventricular assist device (VAD), and clear advantages of the VMS stabilization compared to the SUPG stabilized formulation are presented.

Tables

Table 1. Table 1 : Characteristics of the simulation run for the different discretizations used for the 2D stirrer test case.

method	$n_{\mathrm{GMRES}}$	$n_{\mathrm{NR}}$	${\epsilon }_{det}$	$\max \left|det\left({𝑺}_i\right)-1\right|$
morph-GLS	2808	930	$6.22\cdot {10}^{-5}$	$7.67\cdot {10}^{-6}$
log-morph-SUPG	603	213	$4.23\cdot {10}^{-12}$	$6.03\cdot {10}^{-13}$
log-morph-VMS	596	212	$4.32\cdot {10}^{-12}$	$4.59\cdot {10}^{-13}$

Table 2. Table 2 : Characteristics of the simulation run for the different stabilizations of the log-morph method for the ReinVAD test case.

method	$n_{\mathrm{GMRES}}$	$n_{\mathrm{NR}}$	${\epsilon }_{det}$	$\max \left|det\left({𝑺}_i\right)-1\right|$
log-morph-SUPG	3600	1800	$0.169$	$0.0880$
log-morph-VMS	1180	590	$2.74\cdot {10}^{-8}$	$1.47\cdot {10}^{-9}$

Equations

$$\mathcal{R}\left(\mathbf{S}\right)=\frac{\partial \mathbf{S}}{\partial t}+\left(\mathbf{u}\cdot \nabla \right)\mathbf{S}+\underset{\text{relaxation}}{\underbrace{{\alpha }_1\left(\mathbf{S}-g\left(\mathbf{S}\right)1\right)}}-\underset{\text{elongation}}{\underbrace{{\alpha }_2\left(\mathbf{E}\mathbf{S}+\mathbf{S}\mathbf{E}\right)}}-\underset{\text{rotation}}{\underbrace{{\alpha }_3\left(\mathbf{W}\mathbf{S}-\mathbf{S}\mathbf{W}\right)}}=0.$$

$$D=\frac{L-B}{L+B}.$$

$${\sigma }_{\text{eff}}=\frac{2\mu {\alpha }_{1}D}{\left(1-D^2\right){\alpha }_2},$$

$$\mathcal{R}\left(\mathbf{\psi }\right)=\frac{\partial \mathbf{\psi }}{\partial t}+\left(\mathbf{u}\cdot \nabla \right)\mathbf{\psi }$$

$$\text{with }\mathbf{F}\left(\mathbf{\psi },\mathbf{E}\right)=$$

$$0=\underset{\equiv W\left(\mathbf{\varphi };\mathbf{\psi }\right)}{\underbrace{{\int }_{Q_n}\mathbf{\varphi }:\mathcal{R}\left(\mathbf{\psi }\right)\mathrm{d}Q}}+\underset{\equiv D\left(\mathbf{\varphi };\mathbf{\psi }\right)}{\underbrace{{\int }_{{\Omega }_{n-1}}\mathbf{\varphi }\left(t_{n-1}^{+}\right):\left(\mathbf{\psi }\left(t_{n-1}^{+}\right)-\mathbf{\psi }\left(t_{n-1}^{-}\right)\right)\mathrm{d}\Omega }}$$

$$W\left(\mathbf{\varphi };\mathbf{\psi }\right)=$$

$$\underset{W_n\left(\mathbf{\varphi };\mathbf{\psi }\right)}{\underbrace{-{\int }_{Q_n}\mathbf{\varphi }:\left[{\alpha }_1\frac{3\exp \left(-\mathbf{\psi }\right)}{\mathrm{tr}\left(\exp \left(-\mathbf{\psi }\right)\right)}+{\alpha }_2\mathbf{F}\left(\mathbf{\psi },\mathbf{E}\right)\right]\mathrm{d}Q}}.$$

$$\mathbf{U}={\mathbf{U}}^h+\tilde{\mathbf{U}}.$$

$${\mathcal{S}}_n={\mathcal{S}}_n^h\oplus {\tilde{\mathcal{S}}}_n\text{and}{\mathcal{V}}_n={\mathcal{V}}_n^h\oplus {\tilde{\mathcal{V}}}_n.$$

$${\mathcal{S}}_n^h$$

$${\mathcal{V}}_n^h$$

$$\mathbf{\psi }={\mathbf{\psi }}^h+\tilde{\mathbf{\psi }}$$

$$0$$

$$W_n\left(\bullet ;{\mathbf{\psi }}^h+\tilde{\mathbf{\psi }}\right)\approx W_n\left(\bullet ;{\mathbf{\psi }}^h\right)+{\frac{\partial }{\partial ϵ}{W}_n\left(\bullet ;{\mathbf{\psi }}^h+ϵ\tilde{\mathbf{\psi }}\right)|}_{ϵ=0},$$

$$W_l\left(\tilde{\mathbf{\varphi }};\tilde{\mathbf{\psi }}\right)+{\frac{\partial }{\partial ϵ}{W}_n\left(\tilde{\mathbf{\varphi }};{\mathbf{\psi }}^h+ϵ\tilde{\mathbf{\psi }}\right)|}_{ϵ=0}+D\left(\tilde{\mathbf{\varphi }};\tilde{\mathbf{\psi }}\right)=-W_l\left(\tilde{\mathbf{\varphi }};{\mathbf{\psi }}^h\right)-W_n\left(\tilde{\mathbf{\varphi }};{\mathbf{\psi }}^h\right)-D\left(\tilde{\mathbf{\varphi }};{\mathbf{\psi }}^h\right).$$

$$\tilde{\mathbf{\psi }}=\mathcal{F}\left({\mathbf{\psi }}^h,\mathcal{R}\left({\mathbf{\psi }}^h\right)\right).$$

$$\tilde{\mathbf{\psi }}\approx -\tau \mathcal{R}\left({\mathbf{\psi }}^h\right)\equiv -\tau {\mathcal{R}}^h,$$

$$\tau ={\alpha }_{\tau }{\left({\left(\frac{2}{\Delta t}\right)}^2+\mathbf{u}\cdot \mathbf{G}\mathbf{u}+{\Vert \mathbf{L}\Vert }_2\right)}^{-\frac{1}{2}},$$

$$W_l\left({\mathbf{\varphi }}^h;{\mathbf{\psi }}^h-\tau {\mathcal{R}}^h\right)=W_l\left({\mathbf{\varphi }}^h;{\mathbf{\psi }}^h\right)$$

$$+{\int }_{Q_n}\tau {\alpha }_3{\mathbf{\varphi }}^h:\left[\mathbf{W}{\mathcal{R}}^h-{\mathcal{R}}^h\mathbf{W}\right]\mathrm{d}Q.$$

$$\frac{\partial }{\partial ϵ}\frac{1}{z-ϵ{\mathbf{\psi }}^h}=\frac{1}{z-ϵ{\mathbf{\psi }}^h}{\mathbf{\psi }}^h\frac{1}{z-ϵ{\mathbf{\psi }}^h},$$

$${\frac{\partial }{\partial ϵ}\exp \left(-{\mathbf{\psi }}^h+ϵ\tau {\mathcal{R}}^h\right)|}_{ϵ=0}$$

$$=-\tau \underset{\equiv \mathbf{K}\left({\mathbf{\psi }}^h,{\mathcal{R}}^h\right)}{\underbrace{\frac{1}{2\pi i}{\int }_{\Gamma }{\mathrm{e}}^{-z}\frac{1}{z-{\mathbf{\psi }}^h}{\mathcal{R}}^h\frac{1}{z-{\mathbf{\psi }}^h}\mathrm{d}z}}.$$

$${\frac{\partial }{\partial ϵ}\frac{\exp \left(-{\mathbf{\psi }}^h+ϵ\tau {\mathcal{R}}^h\right)}{\mathrm{tr}\left(\exp \left(-{\mathbf{\psi }}^h+ϵ\tau {\mathcal{R}}^h\right)\right)}|}_{ϵ=0}=-\tau \underset{\equiv {\mathbf{L}}_{{\alpha }_1}\left({\mathbf{\psi }}^h,{\mathcal{R}}^h\right)}{\underbrace{\left(\frac{\mathbf{K}\left({\mathbf{\psi }}^h,{\mathcal{R}}^h\right)}{\mathrm{tr}\left(\exp \left(-{\mathbf{\psi }}^h\right)\right)}-\frac{\exp \left(-{\mathbf{\psi }}^h\right)}{\mathrm{tr}{\left(\exp \left(-{\mathbf{\psi }}^h\right)\right)}^2}\mathrm{tr}\left(\mathbf{K}\left({\mathbf{\psi }}^h,{\mathcal{R}}^h\right)\right)\right)}},$$

$${\frac{\partial }{\partial ϵ}\mathbf{F}\left({\mathbf{\psi }}^h-ϵ\tau {\mathcal{R}}^h,\mathbf{E}\right)|}_{ϵ=0}$$

$$=-\tau \underset{\equiv {\mathbf{L}}_{{\alpha }_2}\left({\mathbf{\psi }}^h,{\mathcal{R}}^h,\mathbf{E}\right)}{\underbrace{\frac{1}{{\left(2\pi i\right)}^2}{\int }_{\Gamma }{\int }_{\Gamma }f\left(z-z^{\prime }\right)\left[\frac{1}{z-{\mathbf{\psi }}^h}{\mathcal{R}}^h\frac{1}{z-{\mathbf{\psi }}^h}\mathbf{E}\frac{1}{z^{\prime }-{\mathbf{\psi }}^h}+\frac{1}{z-{\mathbf{\psi }}^h}\mathbf{E}\frac{1}{z^{\prime }-{\mathbf{\psi }}^h}{\mathcal{R}}^h\frac{1}{z^{\prime }-{\mathbf{\psi }}^h}\right]\mathrm{d}z\mathrm{d}{z}^{\prime }}}.$$

$$0$$

$$+{\int }_{Q_n}\tau \left[\frac{\partial {\mathbf{\varphi }}^h}{\partial t}+\left(\mathbf{u}\cdot \nabla \right){\mathbf{\varphi }}^h\right]:\mathcal{R}\left({\mathbf{\psi }}^h\right)\mathrm{d}Q+{\int }_{Q_n}\tau {\alpha }_3{\mathbf{\varphi }}^h:\left[\mathbf{W}\mathcal{R}\left({\mathbf{\psi }}^h\right)-\mathcal{R}\left({\mathbf{\psi }}^h\right)\mathbf{W}\right]\mathrm{d}Q$$

$$+{\int }_{Q_n}\tau {\mathbf{\varphi }}^h:\left[3{\alpha }_1{\mathbf{L}}_{{\alpha }_1}\left({\mathbf{\psi }}^h,\mathcal{R}\left({\mathbf{\psi }}^h\right)\right)+{\alpha }_2{\mathbf{L}}_{{\alpha }_2}\left({\mathbf{\psi }}^h,\mathcal{R}\left({\mathbf{\psi }}^h\right),\mathbf{E}\right)\right]\mathrm{d}Q.$$

$$\frac{1}{z-{\mathbf{\psi }}^h}=\sum _{i=1}^d\frac{1}{z-{\lambda }_i}{\mathcal{P}}_i.$$

$$\frac{1}{2\pi i}{\int }_{\Gamma }m\left(z\right)\frac{1}{z-\lambda }\mathrm{d}z=m\left(\lambda \right),$$

$$\exp \left(-{\mathbf{\psi }}^h\right)=\frac{1}{2\pi i}{\int }_{\Gamma }{\mathrm{e}}^{-z}\frac{1}{z-{\mathbf{\psi }}^h}\mathrm{d}z=\sum _{i=1}^d{\mathrm{e}}^{-{\lambda }_i}{\mathcal{P}}_i,$$

$$\mathrm{tr}\left(\exp \left(-{\mathbf{\psi }}^h\right)\right)=\sum _{i=1}^d{\mathrm{e}}^{-{\lambda }_i}$$

Full text

Keywords: computational hemodynamics; variational multiscale formulation; log-morphology formulation; finite element method; ventricular assist device

1 Introduction

Computational analysis has become a main step in the development of blood-handling medical devices. Computational Fluid Dynamics (CFD), in particular, helps to reduce the number of expensive and time-consuming experiments during prototype construction. CFD is able to predict the hydraulic performance and the flow behavior within such devices sufficiently and can therefore give useful information towards improvements of their design. This information is very valuable, since the seventh INTERMACS annual report [1] showed that adverse event rates for hemolysis, strokes and renal dysfunction have increased for patients with a ventricular assist device during the last survey period.

The accurate numerical prediction of hemolysis remains a challenging task. An empirical power law model is widely used to estimate the produced plasma-free hemoglobin in medical devices [2, 3]. In such a model, the amount of generated plasma-free hemoglobin is proportional to the shear stress and the duration over which the RBCs are exposed to that stress. The power law parameters are normally determined by a fitting to experimental data obtained by Couette shearing device experiments using human or animal blood samples [4, 5, 6]. This modeling approach, the so-called stress-based hemolysis model, always computes a scalar measure of the shear stress from the flow field, which assumes an instantaneous deformation of the red blood cells.

Since such a stress-based model is not taking the viscoelastic behavior of RBCs into account, Arora *et al. * [7] proposed a strain-based model that estimates the droplet-like deformation of red blood cells in blood flow. The so-called morphology equation is able to account for the relaxation, elongation and rotation of the droplets in the flow. Arora *et al. * [7] used the model in a Lagrangian frame, and Pauli *et al. * [8] proposed its application in an Eulerian frame. Alternative viscoelastic blood damage models are introduced by different authors. Chen and Sharp [9] used a cell threshold model that they calibrated to channel flow experiments to estimate fatal hemolysis. However, such a threshold model cannot be applied to sublethal hemolysis. Arwatz and Smits [10] derived a scalar viscoelastic blood damage model, yet the generalisation to three-dimensional flows is not clear. Other strain-based approaches were proposed by Ezzeldin *et al. * [11] and Sohrabi and Liu [12]. Ezzeldin *et al. * [11] compute the deformation of a high-fidelity RBC model based on a membrane energy composed of four contributions. Sohrabi and Liu [12] use a spring-connected network model for the description of the RBC membrane and use a pore formation model to estimate the generated plasma-free hemoglobin. Both works (Refs. [11, 12]) describe the red blood cell membrane accurately but rely on a Lagrangian description, which can be biased by the choice of the tracers. Furthermore, it is computationally too costly to apply these methods to large-scale simulations such as those required for whole blood pumps.

In this study, we will use the morphology equation to compute the deformation of RBCs. Since the ellipsoidal RBC shape is described by the square root of the eigenvalues of the shape tensor, it is important that these eigenvalues are all positive, i.e., the shape tensor has to be positive-definite in order to describe a physical state. This positive-definiteness can be violated during the simulation, which is also a well-known issue for the related Oldroyd-B model. There, a logarithmic transformation of the conformation tensor fulfills the positive-definiteness by construction, and hence, improves the numerical stability of the method [13, 14]. We will make use of such a transformation for the morphology equation and stabilize the transformed equations with a variational multiscale formulation.

The variational multiscale (VMS) method was introduced by Hughes [15] and Hughes and Stewart [16] as a formulation to derive stabilization terms for the finite element method as analytical corrections from unresolved fine-scale contributions to the governing equations. In addition to the first applications to advection-diffusion and to Helmholtz equations, the VMS concept was applied to various other fields such as, e.g., the turbulence modeling for the incompressible Navier-Stokes equations [17] and also to tensor-based equations such as the Oldroyd-B model [18].

The structure of the paper is as follows: In the next section, we will shortly introduce the morphology equation for the simulation of RBC deformation. We will motivate and present the transformation to a logarithmic shape tensor for enhanced numerical stability and its weak form for a finite-element implementation. The section will be completed with a detailed derivation of the variational multiscale formulation for this log-morph equation. Section 3 will cover the numerical implementation of the newly derived VMS terms as well as the derivation and treatment of their directional derivatives for the usage in a Newton-Raphson solver. Two test cases, a two dimensional stirrer and a state-of-the-art ventricular assist device, will be investigated in Section 4 and the advantages of the proposed log-morph formulation with VMS stabilization terms will be presented.

2 Morphology Equation

The behavior of red blood cells (RBCs) in blood flow is dependent on the shear rates in the fluid. At low shear rates, RBCs tend to form stack-like structures, called rouleaux [19, 20]. A moderate increase of the fluid shear to tens of $\mathrm{s}^{-1}$ breaks these rouleaux and the individual, biconcave RBCs move and tumble in the plasma [21]. A further increase of the shear rate lets the RBCs rotate in the flow and deform to an ellipsoidal shape with a strained membrane [22]; at around $4000\text{\,}\mathrm{s}^{-1}$, pores form in the lipid bilayer to release the membrane stresses, through which hemoglobin is released to the blood plasma [12, 23]. Very high and unphysiological shear rates of about $42\,000\text{\,}\mathrm{s}^{-1}$ can ultimately lead to fatal hemolysis, the complete rupture of the RBC [24]. A recent study showed that the rotation of the RBC’s membrane around the enclosed cytoplasm, the so-called tank treading motion, may not be observed for physiological viscosity ratios [25].

2.1 Droplet Model

Arora et al. [7] proposed in 2004 a droplet model to account for the relaxation, deformation, and rotation of the RBCs, the so-called morphology model. This tensor-based model is able to describe the deformed, ellipsoidal shape of the RBCs. Arora *et al. * used this model in a Lagrangian frame, and included the tank treading motion in their model. Recent advances (cf. Pauli *et al. * [8]) used the tensor-based model in an Eulerian frame without the consideration of the tank-treading motion.

In the morphology model, the ellipsoidal RBC is described by a symmetric, positive-definite $3\times 3$ matrix $\boldsymbol{S}$. The residual $\boldsymbol{\mathcal{R}}$ of the governing equation for the behavior of the RBCs in an external flow field $\boldsymbol{u}$ is given by

[TABLE]

with unit matrix $\boldsymbol{1}$, the strain rate tensor $\boldsymbol{E}=\left(\nabla\boldsymbol{u}+\nabla\boldsymbol{u}^{T}\right)/2$ and the vorticity tensor $\boldsymbol{W}=\left(\nabla\boldsymbol{u}-\nabla\boldsymbol{u}^{T}\right)/2$. The parameter $g\!\left(\boldsymbol{S}\right)=\frac{3I\!I\!I\left(\boldsymbol{S}\right)}{I\!I\left(\boldsymbol{S}\right)}$, with the third invariant $I\!I\!I\!\left(\boldsymbol{S}\right)=\det\!\left(\boldsymbol{S}\right)$ and the second invariant $I\!I\!\left(\boldsymbol{S}\right)=\left(\mathrm{tr}\!\left(\boldsymbol{S}\right)^{2}-\mathrm{tr}\!\left(\boldsymbol{S}^{2}\right)\right)/2$, ensures the conservation of the volume of the RBCs. The parameters $\alpha_{1}=$5\text{,}\mathrm{s}^{-1} and $\alpha_{2}=\alpha_{3}=$4.2298\text{\times}{10}^{-4} were derived by Arora et al. [7] from RBC relaxation and deformation properties.

The square roots of the eigenvalues of $\boldsymbol{S}$ are the semi-axes lengths of the ellipsoid. For a computed shape $\boldsymbol{S}$, one can compute the distortion $D$, a measure of the RBC’s deformation, with the longest and shortest semi-axis, $L$ and $B$ by

[TABLE]

One can show for a simple shear flow that the distortion is a function of the scalar shear stress $\sigma_{\textrm{f}}=2\mu\sqrt{-I\!I\!\left(\boldsymbol{E}\right)}$, where $\mu$ is the blood viscosity. By inverting this relation, we get a formula for an effective shear stress that is acting on the RBC,

[TABLE]

which is generalized and also used for complex flow and transient situations.

2.2 Fully-Implicit Log-Morphology Formulation

The positive definiteness of the morphology tensor $\boldsymbol{S}$ is a necessary condition to describe a physical ellipsoidal shape. Nevertheless, this condition can be numerically violated leading to a diverged simulation. The morphology equation resembles the upper-convected Maxwell model, for which the positive-definiteness of the conformation tensor is also hard to satisfy numerically. Fattal and Kupferman [13] proposed to use a log-conformation tensor — a transformation that satisfies the positive definiteness by design — to circumvent this problem. Knechtges et al. [26] and Knechtges [14] derived a fully-implicit log-conformation formulation, which can be also used with some modifications for the morphology equation.

The governing equation for the log-morphology tensor $\boldsymbol{\psi}=\log\left(\boldsymbol{S}\right)$, or in short, log-morph equation, can be derived as

[TABLE]

representing a double Cauchy integral with the function $f\!\left(z\right)=\frac{z}{\tanh\left(z/2\right)}$ and a suitable contour $\Gamma$, which only encloses the complete spectrum of $\boldsymbol{\psi}$, but not the poles of $f\!\left(z-z^{\prime}\right)$. The derivation of the volume conservation term $g\!\left(\boldsymbol{\psi}\right)=\frac{3}{\mathrm{tr}\left(\exp\left(-\boldsymbol{\psi}\right)\right)}$ is shown in Appendix A.

The arguments of Knechtges [14] hold also for the log-morph equation (4) and a solution of this equation, transformed back, is also a solution of the untransformed morphology equation (1).

Weak Form

We solve the log-morph equation using a space-time finite element method. To this end, we extrude our spatial domain in the time direction and divide the simulation time $t\in\left[0,T\right)$ into subintervals $I_{n}=\left[t_{n-1},t_{n}\right)$. We introduce space-time slabs $Q_{n}$ that are bounded by the spatial domain $\Omega_{n-1}$ at the time-step $t_{n-1}$ and by $\Omega_{n}$ at $t_{n}$. The spatial boundary of the space-time slab is denoted by $P_{n}$. Let the trial solution and weighting function spaces over $Q_{n}$ be given as $\mathcal{S}_{n}$ and $\mathcal{V}_{n}$. The variational form of the log-morph equation on a space-time slab $Q_{n}$ can then be given as: Find $\boldsymbol{\psi}\in\mathcal{S}_{n}$ for the given initial condition $\boldsymbol{\psi}\!\left(t_{0}^{-}\right)=\boldsymbol{\psi}_{0}$ such that $\forall\,\boldsymbol{\phi}\in\mathcal{V}_{n}$

[TABLE]

is satisfied. The second integral $D\!\left(\boldsymbol{\phi};\boldsymbol{\psi}\right)$ is a discontinuous Galerkin term that weakly imposes the continuity of $\boldsymbol{\psi}$ across space-time slabs. The Galerkin term $W\!\left(\boldsymbol{\phi};\boldsymbol{\psi}\right)$ can be decomposed into its linear and nonlinear parts

[TABLE]

2.3 Variational Multiscale Formulation for Log-Morph

In the VMS approach, one assumes that the real solution $\boldsymbol{U}$ to a general problem is composed of a coarse-scale solution $\boldsymbol{U}^{h}$ that can be resolved by the numerical method and an unresolved, fine-scale solution $\widetilde{\boldsymbol{U}}$ (cf. References [15, 17]), i.e.,

[TABLE]

A further assumption is that the trial solution and weighting function spaces can be described as a direct sum of the coarse scale and fine scale spaces (cf. Ref. [17]), i.e.,

[TABLE]

For the log-morph equation, we define the corresponding coarse scale trial solution and weighting function spaces as

[TABLE]

on the finite element mesh $\overline{Q_{n}}$ and with the subset $P_{n}^{D}$ of the spatial space-time slab boundary $P_{n}$ where Dirichlet boundary conditions $\boldsymbol{\psi}=\boldsymbol{g}_{\boldsymbol{\psi}}$ are prescribed. Let us further assume that the fine scale trial functions are $\widetilde{\boldsymbol{\psi}}\in\widetilde{\mathcal{S}}_{n}$ and the corresponding weighting functions are $\widetilde{\boldsymbol{\phi}}\in\widetilde{\mathcal{V}}_{n}$ and that the VMS assumption

[TABLE]

holds. The weak form of the log-morph equation (5) is then decomposed into a set of coupled coarse and fine scale equations

[TABLE]

We can linearize the nonlinear terms in eqs. (12) and (13) using Fréchet derivatives [27]:

[TABLE]

which lets us rearrange eq. (13) to

[TABLE]

This is the weak form of a differential equation for the fine-scale solution $\widetilde{\boldsymbol{\psi}}$, where the left-hand-side only contains operators acting on $\widetilde{\boldsymbol{\psi}}$ and where only terms containing $\boldsymbol{\psi}^{h}$ and $\boldsymbol{\mathcal{R}}^{h}$ occur on the right-hand-side. As in Bazilevs et al. [17], eq. (15) hence tells us that the fine scale solution is a functional of the coarse scale solution and its residual

[TABLE]

Using physical reasoning, one can use a perturbation series for the fine scale solution that is dependent on powers of the coarse scale residual. One can show for the lowest order approximation that

[TABLE]

with the stabilization parameter $\tau$ that we choose according to Shakib et al. [28] as

[TABLE]

where $\boldsymbol{G}_{ij}=\sum_{k}\frac{\partial\xi^{k}}{\partial x^{i}}\frac{\partial\xi^{k}}{\partial x^{j}}$ is the covariant metric tensor mapping to a symmetric111We include the mapping to an equilateral triangle or a regular tetrahedron from Pauli [29] in the definition of the metric tensor. reference element [29] and $\left\|\boldsymbol{L}\right\|_{2}$ is the spectral norm of the linearized source term, which we choose as the Jacobian of the source term (see eq. (34)). We also use a scale factor $\alpha_{\tau}$ for our numerical implementation. In general, Hughes [15] and others [16, 30, 17] showed that the matrix-valued stabilization parameter $\boldsymbol{\tau}$ can be computed with the fine scales Green’s operator.

Inserting the approximation from eq. (17) into eq. (12) closes the system and lets us solve for the coarse scale log-morph tensor $\boldsymbol{\psi}^{h}$. It can be easily deduced by an integration by parts that the VMS formulation includes the SUPG stabilization terms, since

[TABLE]

As mentioned before, the arising nonlinear terms have to be linearized with respect to the coarse scale solution. In order to compute the Fréchet derivatives, we have to use the following relation,

[TABLE]

and we need to compute the derivative of the matrix exponential for the $\alpha_{1}$-term,

[TABLE]

The resulting linearized terms $\boldsymbol{L}_{\alpha_{1}}$ and $\boldsymbol{L}_{\alpha_{2}}$ are:

[TABLE]

and

[TABLE]

It is important to note, that the additional VMS terms are all traceless which ensures the volume conservation (cf. Appendix A).

The resulting VMS stabilized equations are

[TABLE]

3 Numerical Implementation

For the numerical implementation of the log-morph equation with the VMS stabilization, we need to perform some further steps. The discretized equations are going to be solved using a Newton-Raphson algorithm combined with a GMRES solver for the resulting linearized equation system. The Newton-Raphson algorithm requires a further directional derivative with respect to $\delta\boldsymbol{\psi}^{h}$ (also a Fréchet derivative) for the linearization of the system. This derivative has to be computed for every term arising in $W_{l}\!\left(\boldsymbol{\phi}^{h};\boldsymbol{\psi}^{h}-\tau\boldsymbol{\mathcal{R}}^{h}\right)$ and $W_{n}\!\left(\boldsymbol{\phi}^{h};\boldsymbol{\psi}^{h}-\tau\boldsymbol{\mathcal{R}}^{h}\right)$.

In order to evaluate the Cauchy integrals that arise, we use projectors $\boldsymbol{\mathcal{P}}_{i}=\boldsymbol{e}_{i}\boldsymbol{e}_{i}^{T}$ that project onto the one-dimensional subspaces of the eigenvalues $\lambda_{i}$ of $\boldsymbol{\psi}^{h}$ spanned by the corresponding eigenvectors $\boldsymbol{e}_{i}$, i.e.,

[TABLE]

This leads to Cauchy integrals that can be evaluated using the residue theorem, which are similar to

[TABLE]

with an arbitrary function $m\!\left(z\right)$ without a singularity at $z=\lambda$.

3.1 Numerical Evaluation of the VMS Weak Form

The evaluation of the Cauchy integrals is done according to Knechtges [14], whose numerical treatment of the arising prefactors can also be used for our equations. The resulting $\alpha_{1}$-term in the residual can be computed with the matrix exponential

[TABLE]

from which

[TABLE]

immediately follows. The $\alpha_{2}$-term can be evaluated with

[TABLE]

The function $f(x)=\frac{x}{\tanh(x/2)}$ is set to 2 for $x=0$ and can be readily used for values of $x\neq 0$.

Similar considerations yield the new VMS terms in eq. (24), where we have to further evaluate

[TABLE]

and take traces of the form

[TABLE]

with an arbitrary $3\times 3$ matrix $\boldsymbol{A}$ for the determination of the $\boldsymbol{L}_{\alpha_{1}}\!\left(\boldsymbol{\psi}^{h},\boldsymbol{\mathcal{R}}^{h}\right)$ term. We find for the other nonlinear VMS term

[TABLE]

As discussed by Knechtges [14], the prefactor $\frac{\sinh\left(x/2\right)}{(x/2)}$ is used for values $x\neq 0$, and replaced by 1 for $x=0$. The prefactor in the $\boldsymbol{L}_{\alpha_{2}}$ term can be approximated by a Taylor series

[TABLE]

in the vicinity of small denominators, $|x-y|<10^{-2}$. The derivatives are approximated by their Taylor series up to fourth order for small arguments $\left|\nicefrac{{\left(x+y\right)}}{{2}}\right|<10^{-1}$.

3.2 Numerical Evaluation of the Directional Derivatives for the Newton-Raphson Algorithm

As mentioned before, we have to take a further derivative in the direction of $\delta\boldsymbol{\psi}^{h}$ for the linearization in the Newton-Raphson algorithm. For the residual $\boldsymbol{\mathcal{R}}^{h}$ this gives us the Jacobian

[TABLE]

where $\delta\boldsymbol{\psi}^{h}$ has to be used instead of $\tau\boldsymbol{\mathcal{R}}^{h}$ in $\boldsymbol{L}_{\alpha_{1}}$ and $\boldsymbol{L}_{\alpha_{2}}$. This Jacobian is used for the Galerkin part as well as for the linearized VMS terms of the weak form in eq. (24).

For the nonlinear $\alpha_{1}$-term, we have to compute the directional derivative of $\boldsymbol{K}\!\left(\boldsymbol{\psi}^{h},\boldsymbol{\mathcal{R}}^{h}\right)$ in the direction of $\delta\boldsymbol{\psi}^{h}$, which leads to

[TABLE]

The prefactor of the second term containing the residuals is also discussed by Knechtges [14] and Taylor series approximations are used for small denominators $|x-y|<10^{-3}$ and for the arising derivatives for arguments $\left|\nicefrac{{(x+y)}}{{2}}\right|<10^{-3}$. For the complete directional derivative of the $\boldsymbol{L}_{\alpha_{1}}$-term in eq. (22), we have to consider the chain rule.

Accordingly, the directional derivative of $\boldsymbol{L}_{\alpha_{2}}$ becomes

[TABLE]

The first prefactor was already discussed in the previous section, but the second and third prefactors have to be treated with special care. A complete discussion of these prefactors is done in Appendix B.

4 Results

In order to test the newly proposed log-morph equation with VMS stabilization, we investigate two test cases: A simple two-dimensional rotating stirrer in a square box, and a state-of-the-art ventricular assist device (VAD). We consider three different discretizations for our test cases: the untransformed morphology equation with GLS stabilization (morph-GLS, we refer to Ref. [31] for a detailed derivation), the log-morph equation with only SUPG stabilization (log-morph-SUPG), and the log-morph equation with full VMS stabilization (log-morph-VMS). The computations for this section were performed on the supercomputer JURECA at Forschungszentrum Jülich [32].

4.1 Stirrer Test Case

The computational mesh for the stirrer test case consists of $n_{n}=$46,647 nodes and $n_{e}=$92\,262 unstructured triangular elements. A part of the mesh near the beam is shown in Fig. 1 together with the dimensions of the geometry.

The flow solution is computed with the in-house deformable-spatial-domain/stabilized space-time finite element solver XNS. We use a Newtonian blood model with a density of $\rho=$1054\text{,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3} and a viscosity of $\mu=$0.0035\text{\,}\mathrm{P}\mathrm{a}\,\mathrm{s}. We utilize the multiple reference frames (MRF) method (cf. Ref. [33]) to compute a steady-state approximation of the flow field for an angular velocity of $\omega=$50\pi\text{,}\mathrm{s}^{-1}$$. The angular velocity is increased from 0 to the full rotational speed through the first 30 Newton-Raphson iterations $i_{\mathrm{NR}}$ according to

[TABLE]

The steady-state flow field is converged after 34 iterations and is shown in Fig. 2, where we show the streamlines in the corresponding reference frame.

The flow solution is used as an ambient velocity for the computation of the deformation of the RBCs. For this simple test case, we are able to choose a time step size $\Delta t=$0.01\text{,}\mathrm{s} for all three discretizations. We compute a quasi-steady deformation for an exposure time of $1\text{\,}\mathrm{s}$ using a semi-discrete time discretization scheme with a backward differentiation formula of second order (BDF2). Also for the morphology equation, we use an MRF approach. A Newton-Raphson algorithm combined with a preconditioned, restarted GMRES solver is used to solve the nonlinear system. For all three simulations we chose a Krylov space of $10$ and an ILUT-preconditioner with a maximal fill-in of $20$ and a threshold of ${10}^{-4}$. In the morph-GLS case we have to further use an augmented Lagrangian method with a penalty parameter of $\varepsilon_{p}=$10\,000 in order to penalize deviations from the initial droplet volume (cf. Ref. [31]).

A comparison of the instantaneous shear stress $\sigma_{\mathrm{f}}$ and the effective shear stress $\sigma_{\mathrm{eff}}$ is shown in Fig. 3(a) and Fig. 3(b).

It can be noted that the morphology equation predicts a one order of magnitude smaller shear stress acting on the RBCs compared to the instantaneous fluid stresses. As can be seen in Figure 3(c), there is hardly any difference between the morph-GLS and the log-morph-SUPG result. The effective shear stress of the log-morph-VMS simulation shows only minor differences in the inner stirrer region, where the stresses are two orders of magnitude smaller compared to the maximum value. The advantage and the enhanced numerical stability of the logarithmic transformation can be seen in Table 1, where we observe a very significant decrease in both the GMRES and Newton-Raphson iterations needed to reach a converged solution.

We use the deviation of the determinant of the shape tensor from 1 over the whole domain, i.e.,

[TABLE]

as a measure for the quality of the volume conservation. Also, this measure illustrates the superiority of the log-morph formulation by a dramatic improvement of the volume conservation, as can be seen in Table 1. For this simple test case, we cannot see any considerable differences in the log-morph-SUPG and the log-morph-VMS formulation. As a final remark, it should be mentioned that the stability of the log-morph formulation allows us to increase the time step size to values of the order of $0.1\text{\,}\mathrm{s}$ without any problems with convergence or the volume conservation.

4.2 State-of-the-Art VAD

The second test case we investigate is a state-of-the-art VAD in preclinical testing developed by the ReinVAD GmbH in Aachen. The computational mesh consists of $n_{n}=$4.74\text{,}\mathrm{M} nodes and $n_{e}=$27.4\text{\,}\mathrm{M} unstructured tetrahedral elements. We use a boundary layer mesh near the no-slip walls of a total thickness of $500\text{\,}\mu\mathrm{m}$ and seven layers with a growth rate of $1.2$, which is compressed in regions with small gap widths. We introduce an MRF interface that encloses the impeller and lies completely in the fluid volume. A part of the computational mesh together with the interface is shown in Fig. 4.

We analyze the pump for an impeller angular velocity of $2400\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$ and a flow rate of $5\text{\,}\mathrm{L}\mathrm{/}\mathrm{m}\mathrm{i}\mathrm{n}$ at the inflow, using again the Newtonian blood model. We compute the steady blood flow with the commercial flow solver Altair AcuSolve using the MRF method and the SST k-$\omega$ turbulence model. The flow solution is depicted for a slice in the middle of the impeller blades in Fig. 5.

The flow field is again used as an ambient velocity for the morphology estimation. For this complex geometry, the untransformed morphology equation, i.e., the morph-GLS method, is not able to give a converged solution. Even with the augmented Lagrangian method, a discontinuity capturing technique, and a time step size as small as $\Delta t=${10}^{-4}\text{,}\mathrm{s}$$, the volume conservation cannot be satisfied and negative eigenvalues in the shape tensor $\boldsymbol{S}$ occur. Nevertheless, the log-morph simulations are able to give results using a scale factor of $\alpha_{\tau}=2$ for the stabilization parameter $\tau$ (eq. (18)) and utilizing a discontinuity capturing similar to the proposition of Shakib *et al. * [28]. The additional discontinuity capturing term in eq. (24) becomes

[TABLE]

with another scale factor $\alpha_{\mathrm{DC}}$, the contravariant metric tensor $\boldsymbol{G}^{-1}$, and a numerical diffusion defined as

[TABLE]

We choose a discontinuity capturing scale factor of $\alpha_{\mathrm{DC}}=$0.05 to ensure convergence. Furthermore, we use the MRF method and a Krylov space of $50$, an ILUT fill-in of $75$ with a threshold of ${10}^{-4}$, and a time step size of $\Delta t=${10}^{-2}\text{\,}\mathrm{s} and simulate for $1.5\text{\,}\mathrm{s}$ physical time. At the inflow, we prescribe a boundary condition for a fully developed pipe flow in axial direction (cf. Ref. [31]).

Figure 6(a) shows the instantaneous shear stress $\sigma_{\mathrm{f}}$ computed from the flow solution on a slice in the middle of the impeller region.

If compared with the effective shear stress $\sigma_{\mathrm{eff}}$ in Fig. 6(b), which is computed from the RBC’s shape estimated by the log-morph-VMS simulation, we see a one order of magnitude lower stress. The instantaneous shear stress is especially high at the impeller blades and their tips in contrast to the effective shear stress. This behavior can be explained with the short exposure times in these regions, leading to relatively low effective stresses acting on the RBCs.

In the line plot in Fig. 6(c) the log-morph-SUPG and the log-morph-VMS method are compared. Both methods predict very similar results in most of the domain. However, close to the no-slip walls at the rotor blades (at a distance of $0.88\text{\,}\mathrm{c}\mathrm{m}$ and $1.32\text{\,}\mathrm{c}\mathrm{m}$) the log-morph-VMS formulation predicts significantly smaller shear stresses compared to the log-morph-SUPG method.

During our simulations, we found convergence issues for the log-morph-SUPG method. The residual fluctuated around a value of $3.7\text{\times}{10}^{-10}$ and hence, did not reach the residual threshold set to ${10}^{-10}$, always using the maximum number of Newton-Raphson iterations of 12 per time step. Another indication for convergence issues is the rather poor volume conservation shown in Table 2.

In contrast, the log-morph-VMS method shows a good convergence and volume conservation using the same solver parameters.

Another interesting comparison is the distribution of the instantaneous shear stress $\sigma_{\mathrm{f}}$ and the effective shear stress $\sigma_{\mathrm{eff}}$ for the two different discretizations on the impeller. Fig. 7 depicts the distribution on top and bottom of the impeller.

It can be noted that the distribution of instantaneous and effective shear stress shows similarities, such as the distribution pattern on the impeller table or in the channels of the hydrodynamic bearing. However, also clear differences are observed, especially at the impeller tips, where the peak values of $\sigma_{\mathrm{f}}$ are not present for $\sigma_{\mathrm{eff}}$. In general, the morphology equation predicts stress values about two orders of magnitude smaller than the instantaneous fluid stresses.

For the discontinuity capturing, we can only use an $\alpha_{\mathrm{DC}}=0.05$ to obtain convergence for the transient log-morph simulations, although we believe that a higher scale factor $\alpha_{\mathrm{DC}}$ is necessary to reduce all oscillations at sharp inter-element discontinuities. However, for the converged log-morph-VMS solution, we are able to compute a restart solution for one more time step with a discontinuity capturing scale factor of $\alpha_{\mathrm{DC}}=$1.0$$, which adds a little bit more numerical diffusion and further decreases peak values. For the log-morph-SUPG simulation, this restart technique did not lead to convergence. From Fig. 7(c), it can be observed that log-morph-VMS results in a smoother effective shear stress distribution with smaller peak values compared to the log-morph-SUPG case. Thus, the VMS formulation helps to reduce oscillations near the impeller walls. Although we present the effective shear stress distribution for the restarted log-morph-VMS method, the log-morph-VMS result before the restart shows already fewer oscillations than the log-morph-SUPG method.

5 Discussion

For the simple stirrer test case, the numerical behavior of the log-morph-SUPG and the log-morph-VMS method show a comparable performance. For the complex VAD test case, though, the VMS-stabilized log-morphology formulation reveals superior convergence behavior and superior stabilization properties at the rotating impeller compared to the SUPG variant. Evidently, this behavior is due to the missing source term stabilization of the SUPG method. The rapidly decreasing velocity magnitudes in the thin boundary layer elements close to the no-slip impeller walls lead to dominating source term contributions. Hence, the VMS-stabilized formulation that by construction includes a source term stabilization is advantageous.

Although the presented VMS method leads to a significant improvement, further developments would be desirable. We expect that an enhanced framework for discontinuity capturing could further improve the convergence rate. Together with an investigation of the definition of the stabilization parameter, this could lead to a better formulation that does not need a scaling of the stabilization terms for complex geometries. Another interesting topic would be the investigation of higher order approximations for the highly nonlinear VMS terms; here, we did only consider the first Fréchet derivatives in the linearization.

6 Conclusion

The aim of this paper was to present the application of the variational multiscale formalism to a tensor-based red blood cell deformation model. We used a logarithmic transformation of the shape tensor in the underlying morphology equation to enhance the numerical stability. This introduces highly nonlinear terms in the resulting log-morph equation. The VMS method is used as a general framework for stabilization of (nonlinear) partial differential equations. Its resulting terms for the log-morph equation are derived with a Fréchet derivative approach for the linearization of the nonlinear terms. To our best knowledge, this is the first application of the VMS formulation to such a highly nonlinear tensor model.

For a simple two dimensional stirrer test case, we found good agreement between the untransformed morphology equation and the log-morph equation with SUPG and VMS stabilization. The enhanced numerical stability of the log-morph equation is clearly observed for this simple test case. Furthermore, we successfully applied the log-morph-VMS method to a state-of-the-art ventricular assist device for which the untransformed morphology equation is not able to converge. Our studies showed that the VMS stabilization method leads to an improved numerical stability for complex test cases compared to an SUPG stabilization scheme.

Acknowledgments

We gratefully acknowledge the fruitful discussions with Philipp Knechtges, without which this paper would not have been possible. We would also like to thank the ReinVAD GmbH for providing us with the investigated VAD geometry data. The authors gratefully acknowledge the computing time granted through JARA-HPC on the supercomputer JURECA at Forschungszentrum Jülich.

Appendix A Derivation of the Volume Conservation Term

The volume conservation for the morphology equation can be derived using the determinant of $\boldsymbol{S}$, i.e., $\det\!\left(\boldsymbol{S}\right)$, which is proportional to the volume of the described ellipsoid. This translates to

[TABLE]

for the logarithmic shape tensor $\boldsymbol{\psi}$. To preserve the volume means mathematically that

[TABLE]

This can be used in eq. (4) and one can show for an incompressible fluid $\left(\mathrm{tr}\!\left(\boldsymbol{E}\right)=0\right)$ that

[TABLE]

and that

[TABLE]

which leads to

[TABLE]

and finally results in

[TABLE]

with the number of space dimensions $d$.

Appendix B Numerical Treatment of the VMS Prefactors

The second prefactor in the directional derivative of the $\boldsymbol{L}_{\alpha_{2}}$-term in eq. (36) is numerically not accurate in the vicinity of small denominators. It can be rewritten with difference quotients using $x=\lambda_{i}-\lambda_{l}$, $y=\lambda_{j}-\lambda_{l}$, and $z=\lambda_{k}-\lambda_{l}$ as

[TABLE]

where the difference quotients are approximated by the Taylor series (33) for values $|x-y|<10^{-2}$ and the arising derivatives are replaced by their Taylor series up to fourth order for $\left|\nicefrac{{(x+y)}}{{2}}\right|<10^{-1}$.

The third prefactor can be approximated by applying the Taylor expansion twice when both denominators are small, yielding

[TABLE]

which is again used when $\left|\lambda_{i}-\lambda_{j}\right|<10^{-2}$ and where the derivatives are approximated for $\left|\nicefrac{{\left(\lambda_{k}+\lambda_{l}-\lambda_{i}-\lambda_{j}\right)}}{{2}}\right|<10^{-1}$.