Ambit – A FEniCS-based cardiovascular multi-physics solver

Author:

Dr.-Ing. Marc Hirschvogel

Preface

Ambit [Hir24] is an open-source multi-physics finite element solver written in Python, supporting solid and fluid mechanics, fluid-structure interaction (FSI), and lumped-parameter models. It is tailored towards solving problems in cardiac mechanics, but may also be used for more general nonlinear finite element analysis. The code encompasses re-implementations and generalizations of methods developed by the author for his PhD thesis [Hir19] and beyond. Ambit makes use of the open-source finite element library FEniCS/dolfinx (https://fenicsproject.org) [LMW12] along with the linear algebra package PETSc (https://petsc.org) [BAA+22], hence guaranteeing a state-of-the-art finite element and linear algebra backend. It is constantly updated to ensure compatibility with a recent dolfinx development version. I/O routines are designed such that the user only needs to provide input files that define parameters through Python dictionaries, hence no programming or in-depth knowledge of any library-specific syntax is required.

Ambit provides general nonlinear (compressible or incompressible) finite strain solid dynamics [Hol00], implementing a range of hyperelastic, viscous, and active material models. Specifically, the well-known anisotropic Holzapfel-Ogden [HO09] and Guccione models [GCM95] for structural description of the myocardium are provided, along with a bunch of other models. It further implements strain- and stress-mediated volumetric growth models [GoktepeAPK10] that allow to model (maladaptive) ventricular shape and size changes. Inverse mechanics approaches to imprint loads into a reference state are implemented using the so-called prestressing method [GForsterW10] in displacement formulation [SG21].

Furthermore, fluid dynamics in terms of incompressible Navier-Stokes/Stokes equations – either in Eulerian or Arbitrary Lagrangian-Eulerian (ALE) reference frames – are implemented. Taylor-Hood elements or equal-order approximations with SUPG/PSPG stabilization [TO00] can be used. Multiphase fluid systems governed by the coupled Navier-Stokes Cahn-Hilliard equations in mass-averaged velocity formulation are also supported [BtE26],:cite:p:ten-eikelder2024, with straightforward extension to ALE formulations.

A variety of reduced 0D lumped models targeted at blood circulation modeling are implemented, including 3- and 4-element Windkessel models [WLW09] as well as closed-loop full circulation [HBJ+17] and coronary flow models [ALA+16].

Monolithic fluid-solid interaction (FSI) [NMK+11] in ALE formulation using a Lagrange multiplier field is supported, along with coupling of 3D and 0D models (solid or fluid with 0D lumped circulation systems) such that cardiovascular simulations with realistic boundary conditions can be performed.

Implementations for a recently proposed novel physics- and projection-based model reduction for FSI, denoted as fluid-reduced-solid interaction (FrSI) [HBBN24], are provided, along with POD-based Galerkin model reduction techniques [FACC14] using full or boundary subspaces.

The nonlinear (single- or multi-field) problems are solved with a customized Newton solver with PTC [GKL09] adaptivity in case of divergence, providing robustness for numerically challenging problems. Linear solvers and preconditioners can be chosen from the PETSc repertoire, and specific block preconditioners are made available for coupled problems.

Avenues for future functionality include cardiac electrophysiology, scalar transport, or finite strain plasticity.

In the following, a brief description of the supported problem types is given, including the strong and weak form of the underlying equations as well as the discrete assembled systems that are solved.

Examples of input files for the respective problem types can be found in the folder demos (with detailed setup descriptions) or amogst the test cases in the folder tests.

This documentation is structured as follows. In sec. 2, instructions on how to install and use Ambit are given. The relevant supported physics models are described in sec. 4. Demos are presented in sec. 5.

Installation

In order to use Ambit, you need to install FEniCSx (https://github.com/FEniCS/dolfinx#installation) (latest Ambit-compatible tested dolfinx release version: v0.10.0.post5).

Ambit can then be installed using pip, either the current release

python3 -m pip install ambit-fe

or latest development version:

python3 -m pip install git+https://github.com/marchirschvogel/ambit.git

Alternatively, you can pull a pre-built Docker image with FEniCSx and Ambit installed:

docker pull ghcr.io/marchirschvogel/ambit:latest

If a Docker image for development is desired, the following image contains all dependencies needed to install and run Ambit:

docker pull ghcr.io/marchirschvogel/ambit:devenv

Ambit input

Here, a minimal Ambit input file is shown, exemplarily for a single-field problem. The mandatory parameter dictionaries to provide are input parameters (IO), global control parameters (CTR), time-integration parameters (TME), solver parameters (SOL), finite element parameters (FEM), and constitutive/material parameters (MAT). For multi-physics problems, each field needs individual time, finite element, and constitutive parameters.

#!/usr/bin/env python3

# Minimal input file for an elastodynamics problem

import ambit_fe
import numpy as np

def main():

    # Input/output
    IO = {"problem_type"        : "solid",                   # type of physics to solve
          "mesh_domain"         : "/path/mesh_d.xdmf",       # path to domain mesh
          "mesh_boundary"       : "/path/mesh_b.xdmf",       # ath to boundary mesh
          "meshfile_type"       : "HDF5",                    # encoding (HDF5 or ASCII)
          "write_results_every" : 1,                         # step frequency for output
          "output_path"         : "/path/...",               # path to output to
          "results_to_write"    : ["displacement",
                                   "vonmises_cauchystress"], # results to output
          "simname"             : "my_results_name"}         # midfix of output name

    # Global control
    CTR = {"maxtime"            : 1.0,  # maximum simulation time
           "dt"                 : 0.01} # time step size

    # Time discretization
    TME = {"timint"             : "genalpha", # time integration: Generalized-alpha
           "rho_inf_genalpha"   : 1.0}        # spectral radius of Gen-alpha scheme

    # Solver
    SOL = {"solve_type"         : "direct", # direct linear solver
           "tol_res"            : 1.0e-8,   # residual tolerance
           "tol_inc"            : 1.0e-8}   # increment tolerance

    # Finite element discretization
    FEM = {"order_disp"         : 1, # FEM degree for displacement field
           "quad_degree"        : 2} # quadrature scheme degree

    # Time curves
    class TC:
        # user defined load curves, to be used in boundary conditions (BC)
        def tc1(self, t):
            load_max = 5.0
            return load_max * np.sin(t)

    # Materials
    MAT = {"MAT1" : {"neohooke_dev" : {"mu" : 10.0},      # isochoric NeoHookean material
                     "ogden_vol"    : {"kappa" : 1.0e3},  # volumetric Ogden material
                     "inertia"      : {"rho0" : 1.0e-6}}} # density

    # Boundary conditions
    BC = {"dirichlet" : [{"id" : [<SURF_IDs>], # list of surfaces for Dirichlet BC
                          "dir" : "all",       # all directions
                          "val" : 0.0}],       # set to zero
          "neumann"   : [{"id" : [<SURF_IDs>], # list of surfaces for Neumann BC
                          "dir" : "xyz_ref",   # in cartesian reference directions
                          "curve" : [1,0,0]}]} # load in x-direction controlled by curve #1 (see time curves)

    # Problem setup
    problem = ambit_fe.ambit_main.Ambit(io_params=IO, ctrl_params=CTR, time_params=TME, solver_params=SOL, fem_params=FEM, constitutive_params=MAT, boundary_conditions=BC, time_curves=TC)

    # Run: solve the problem
    problem.solve_problem()

if __name__ == "__main__":

    main()

Physics Models

Solid mechanics

– Example: Sec. 5.1 and demos/solid

– problem_type : "solid"

– Solid mechanics are formulated in a Total Lagrangian frame

Strong form

Displacement-based

– Primary variable: displacement \(\boldsymbol{u}\)

(1)\[\begin{split}\begin{aligned} \nabla_{0} \cdot \boldsymbol{P}(\boldsymbol{u},\boldsymbol{v}(\boldsymbol{u})) + \hat{\boldsymbol{b}}_{0} &= \rho_{0} \boldsymbol{a}(\boldsymbol{u}) &&\text{in} \; \mathit{\Omega}_{0} \times [0, T], \\ \boldsymbol{u} &= \hat{\boldsymbol{u}} &&\text{on} \; \mathit{\Gamma}_{0}^{\mathrm{D}} \times [0, T],\\ \boldsymbol{t}_{0} = \boldsymbol{P}\boldsymbol{n}_{0} &= \hat{\boldsymbol{t}}_{0} &&\text{on} \; \mathit{\Gamma}_{0}^{\mathrm{N}} \times [0, T],\\ \boldsymbol{u}(\boldsymbol{x}_{0},0) &= \hat{\boldsymbol{u}}_{0}(\boldsymbol{x}_{0}) &&\text{in} \; \mathit{\Omega}_{0},\\ \boldsymbol{v}(\boldsymbol{x}_{0},0) &= \hat{\boldsymbol{v}}_{0}(\boldsymbol{x}_{0}) &&\text{in} \; \mathit{\Omega}_{0}, \end{aligned}\end{split}\]

Incompressible mechanics

– Primary variables: displacement \(\boldsymbol{u}\) and pressure \(p\)

(2)\[\begin{split}\begin{aligned} \nabla_{0} \cdot \boldsymbol{P}(\boldsymbol{u},p,\boldsymbol{v}(\boldsymbol{u})) + \hat{\boldsymbol{b}}_{0} &= \rho_{0} \boldsymbol{a}(\boldsymbol{u}) &&\text{in} \; \mathit{\Omega}_{0} \times [0, T], \\ J(\boldsymbol{u})-1 &= 0 &&\text{in} \; \mathit{\Omega}_{0} \times [0, T], \\ \boldsymbol{u} &= \hat{\boldsymbol{u}} &&\text{on} \; \mathit{\Gamma}_{0}^{\mathrm{D}} \times [0, T],\\ \boldsymbol{t}_{0} = \boldsymbol{P}\boldsymbol{n}_{0} &= \hat{\boldsymbol{t}}_{0} &&\text{on} \; \mathit{\mathit{\Gamma}}_{0}^{\mathrm{N}} \times [0, T],\\ \boldsymbol{u}(\boldsymbol{x}_{0},0) &= \hat{\boldsymbol{u}}_{0}(\boldsymbol{x}_{0}) &&\text{in} \; \mathit{\mathit{\Omega}}_{0},\\ \boldsymbol{v}(\boldsymbol{x}_{0},0) &= \hat{\boldsymbol{v}}_{0}(\boldsymbol{x}_{0}) &&\text{in} \; \mathit{\mathit{\Omega}}_{0}, \end{aligned}\end{split}\]

Weak form

Displacement-based

Find \(\boldsymbol{u}\in\boldsymbol{\mathcal{V}}_{u}^{h,D}\) such that

(3)\[\begin{aligned} r(\boldsymbol{u};\delta\boldsymbol{u}) := \delta \mathcal{W}_{\mathrm{kin}}(\boldsymbol{u};\delta\boldsymbol{u}) + \delta \mathcal{W}_{\mathrm{int}}(\boldsymbol{u};\delta\boldsymbol{u}) - \delta \mathcal{W}_{\mathrm{ext}}(\boldsymbol{u};\delta\boldsymbol{u}) = 0, \end{aligned}\]

for all \(\delta\boldsymbol{u}\in\boldsymbol{\mathcal{V}}_{u}^{h}\)

– Kinetic virtual work:

(4)\[\begin{aligned} \delta \mathcal{W}_{\mathrm{kin}}(\boldsymbol{u};\delta\boldsymbol{u}) &= \int\limits_{\mathit{\Omega}_{0}} \rho_{0}\,\boldsymbol{a}(\boldsymbol{u}) \cdot \delta\boldsymbol{u} \,\mathrm{d}V \end{aligned}\]

– Internal virtual work:

(5)\[\begin{aligned} \delta \mathcal{W}_{\mathrm{int}}(\boldsymbol{u};\delta\boldsymbol{u}) &= \int\limits_{\mathit{\Omega}_{0}} \boldsymbol{P}(\boldsymbol{u},\boldsymbol{v}(\boldsymbol{u})) : \nabla_{0} \delta\boldsymbol{u} \,\mathrm{d}V = \int\limits_{\mathit{\Omega}_{0}} \boldsymbol{S}(\boldsymbol{u},\boldsymbol{v}(\boldsymbol{u})) : \frac{1}{2}\delta\boldsymbol{C}(\boldsymbol{u}) \,\mathrm{d}V \end{aligned}\]

– External virtual work:

• conservative Neumann load:

  (6)\[\begin{aligned} \delta \mathcal{W}_{\mathrm{ext}}(\delta\boldsymbol{u}) &= \int\limits_{\mathit{\Gamma}_{0}^{\mathrm{N}}} \hat{\boldsymbol{t}}_{0}(t) \cdot \delta\boldsymbol{u} \,\mathrm{d}A \end{aligned}\]

• Neumann pressure load in current normal direction:

  (7)\[\begin{aligned} \delta \mathcal{W}_{\mathrm{ext}}(\boldsymbol{u};\delta\boldsymbol{u}) &= -\int\limits_{\mathit{\Gamma}_{0}^{\mathrm{N}}} \hat{p}(t)\,J \boldsymbol{F}^{-\mathrm{T}}\boldsymbol{n}_{0} \cdot \delta\boldsymbol{u} \,\mathrm{d}A \end{aligned}\]

• general Neumann load in current direction:

  (8)\[\begin{aligned} \delta \mathcal{W}_{\mathrm{ext}}(\boldsymbol{u};\delta\boldsymbol{u}) &= \int\limits_{\mathit{\Gamma}_0} J\boldsymbol{F}^{-\mathrm{T}}\,\hat{\boldsymbol{t}}_{0}(t) \cdot \delta\boldsymbol{u} \,\mathrm{d}A \end{aligned}\]

• body force:

  (9)\[\begin{aligned} \delta \mathcal{W}_{\mathrm{ext}}(\delta\boldsymbol{u}) &= \int\limits_{\mathit{\Omega}_{0}} \hat{\boldsymbol{b}}_{0}(t) \cdot \delta\boldsymbol{u} \,\mathrm{d}V \end{aligned}\]

• generalized Robin condition:

  (10)\[\begin{aligned} \delta \mathcal{W}_{\mathrm{ext}}(\boldsymbol{u};\delta\boldsymbol{u}) &= -\int\limits_{\mathit{\Gamma}_{0}^{\mathrm{R}}} \left[k\,\boldsymbol{u} + c\,\boldsymbol{v}(\boldsymbol{u})\right] \cdot \delta\boldsymbol{u}\,\mathrm{d}A \end{aligned}\]

• generalized Robin condition in reference surface normal direction:

  (11)\[\begin{aligned} \delta \mathcal{W}_{\mathrm{ext}}(\boldsymbol{u};\delta\boldsymbol{u}) &= -\int\limits_{\mathit{\Gamma}_{0}^{\mathrm{R}}} (\boldsymbol{n}_0 \otimes \boldsymbol{n}_0)\left[k\,\boldsymbol{u} + c\,\boldsymbol{v}(\boldsymbol{u})\right] \cdot \delta\boldsymbol{u}\,\mathrm{d}A \end{aligned}\]

Incompressible mechanics: 2-field displacement and pressure variables

Find \((\boldsymbol{u},p) \in \boldsymbol{\mathcal{V}}_{u}^{h,D} \times \mathcal{V}_{p}^{h,D}\) such that

(12)\[\begin{split}\begin{aligned} r_u(\boldsymbol{u},p;\delta\boldsymbol{u}) &:= \delta \mathcal{W}_{\mathrm{kin}}(\boldsymbol{u};\delta\boldsymbol{u}) + \delta \mathcal{W}_{\mathrm{int}}(\boldsymbol{u},p;\delta\boldsymbol{u}) - \delta \mathcal{W}_{\mathrm{ext}}(\boldsymbol{u};\delta\boldsymbol{u}) = 0, \\ r_p(\boldsymbol{u};\delta p) &:= \delta \mathcal{W}_{\mathrm{pres}}(\boldsymbol{u};\delta p) = 0, \end{aligned}\end{split}\]

for all \((\delta\boldsymbol{u},\delta p) \in \boldsymbol{\mathcal{V}}_{u}^{h} \times \mathcal{V}_{p}^{h}\)

– Kinetic virtual work: ([equation-deltaw-kin])

– Internal virtual work:

(13)\[\begin{aligned} \delta \mathcal{W}_{\mathrm{int}}(\boldsymbol{u},p;\delta\boldsymbol{u}) &= \int\limits_{\mathit{\Omega}_{0}} \boldsymbol{P}(\boldsymbol{u},p,\boldsymbol{v}(\boldsymbol{u})) : \nabla_{0} \delta\boldsymbol{u} \,\mathrm{d}V = \int\limits_{\mathit{\Omega}_{0}} \boldsymbol{S}(\boldsymbol{u},p,\boldsymbol{v}(\boldsymbol{u})) : \frac{1}{2}\delta\boldsymbol{C}(\boldsymbol{u}) \,\mathrm{d}V \end{aligned}\]

– Pressure virtual work:

(14)\[\begin{aligned} \delta \mathcal{W}_{\mathrm{pres}}(\boldsymbol{u};\delta p) &= \int\limits_{\mathit{\Omega}_{0}} (J(\boldsymbol{u}) - 1) \,\delta p \,\mathrm{d}V \end{aligned}\]

Material models

– hyperelastic material models

\[\begin{aligned} \boldsymbol{S} = 2\frac{\partial\mathit{\Psi}}{\partial \boldsymbol{C}} \end{aligned}\]

- MAT "neohooke_dev"

\[\begin{aligned} \mathit{\Psi} &= \frac{\mu}{2}\left(\bar{I}_C - 3\right) \end{aligned}\]

- MAT "holzapfelogden_dev"

\[\begin{split}\begin{aligned} \mathit{\Psi} &= \frac{a_0}{2b_0}\left(e^{b_0(\bar{I}_C - 3)} - 1\right) + \sum\limits_{i\in\{f,s\}}\frac{a_i}{2b_i}\left(e^{b_i(I_{4,i}-1)^2}-1\right) + \frac{a_{fs}}{2b_{fs}}\left(e^{b_{fs}I_{8}^2} - 1\right), \\ & I_{4,f} = \boldsymbol{f}_0 \cdot \boldsymbol{C}\boldsymbol{f}_0, \quad I_{4,s} = \boldsymbol{s}_0 \cdot \boldsymbol{C}\boldsymbol{s}_0, \quad I_8 = \boldsymbol{f}_0 \cdot \boldsymbol{C}\boldsymbol{s}_0 \end{aligned}\end{split}\]

– viscous material models

\[\begin{aligned} \boldsymbol{S} = 2\frac{\partial\mathit{\Psi}_{\mathrm{v}}}{\partial \dot{\boldsymbol{C}}} \end{aligned}\]

Time integration

– time scheme timint : "static"

\[\begin{aligned} \delta \mathcal{W}_{\mathrm{int}}(\boldsymbol{u}_{n+1};\delta\boldsymbol{u}) - \delta \mathcal{W}_{\mathrm{ext}}(\boldsymbol{u}_{n+1};\delta\boldsymbol{u}) = 0, \quad \forall \; \delta\boldsymbol{u}\end{aligned}\]

– Generalized-alpha time scheme timint : "genalpha"

\[\begin{split}\begin{aligned} \boldsymbol{v}_{n+1} &= \frac{\gamma}{\beta\Delta t}(\boldsymbol{u}_{n+1}-\boldsymbol{u}_{n}) - \frac{\gamma-\beta}{\beta} \boldsymbol{v}_{n} - \frac{\gamma-2\beta}{2\beta}\Delta t\,\boldsymbol{a}_{n} \\ \boldsymbol{a}_{n+1} &= \frac{1}{\beta\Delta t^2}(\boldsymbol{u}_{n+1}-\boldsymbol{u}_{n}) - \frac{1}{\beta\Delta t} \boldsymbol{v}_{n} - \frac{1-2\beta}{2\beta}\boldsymbol{a}_{n} \end{aligned}\end{split}\]

• option eval_nonlin_terms : "midpoint":

(15)\[\begin{split}\begin{aligned} \boldsymbol{u}_{n+1-\alpha_{\mathrm{f}}} &= (1-\alpha_{\mathrm{f}})\boldsymbol{u}_{n+1} + \alpha_{\mathrm{f}} \boldsymbol{v}_{n} \\ \boldsymbol{v}_{n+1-\alpha_{\mathrm{f}}} &= (1-\alpha_{\mathrm{f}})\boldsymbol{v}_{n+1} + \alpha_{\mathrm{f}} \boldsymbol{v}_{n} \\ \boldsymbol{a}_{n+1-\alpha_{\mathrm{m}}} &= (1-\alpha_{\mathrm{m}})\boldsymbol{a}_{n+1} + \alpha_{\mathrm{m}} \boldsymbol{a}_{n} \end{aligned}\end{split}\]

\[\begin{aligned} \delta \mathcal{W}_{\mathrm{kin}}(\boldsymbol{a}_{n+1-\alpha_{m}};\delta\boldsymbol{u}) + \delta \mathcal{W}_{\mathrm{int}}(\boldsymbol{u}_{n+1-\alpha_{f}};\delta\boldsymbol{u}) - \delta \mathcal{W}_{\mathrm{ext}}(\boldsymbol{u}_{n+1-\alpha_{f}};\delta\boldsymbol{u}) = 0, \quad \forall \; \delta\boldsymbol{u}\end{aligned}\]

• option eval_nonlin_terms : "trapezoidal":

\[\begin{split}\begin{aligned} &(1-\alpha_{\mathrm{m}})\,\delta \mathcal{W}_{\mathrm{kin}}(\boldsymbol{a}_{n+1};\delta\boldsymbol{u}) + \alpha_{\mathrm{m}}\,\delta \mathcal{W}_{\mathrm{kin}}(\boldsymbol{a}_{n};\delta\boldsymbol{u}) + \\ & (1-\alpha_{\mathrm{f}})\,\delta \mathcal{W}_{\mathrm{int}}(\boldsymbol{u}_{n+1};\delta\boldsymbol{u}) + \alpha_{\mathrm{f}}\,\delta \mathcal{W}_{\mathrm{int}}(\boldsymbol{u}_{n};\delta\boldsymbol{u}) - \\ & (1-\alpha_{f})\,\delta \mathcal{W}_{\mathrm{ext}}(\boldsymbol{u}_{n+1};\delta\boldsymbol{u}) - \alpha_{\mathrm{f}}\,\delta \mathcal{W}_{\mathrm{ext}}(\boldsymbol{u}_{n};\delta\boldsymbol{u}) = 0, \quad \forall \; \delta\boldsymbol{u}\end{aligned}\end{split}\]

– One-Step-theta time scheme timint : "ost"

\[\begin{split}\begin{aligned} \boldsymbol{v}_{n+1} &= \frac{1}{\theta\Delta t}(\boldsymbol{u}_{n+1}-\boldsymbol{u}_{n}) - \frac{1-\theta}{\theta} \boldsymbol{v}_{n} \\ \boldsymbol{a}_{n+1} &= \frac{1}{\theta^2\Delta t^2}(\boldsymbol{u}_{n+1}-\boldsymbol{u}_{n}) - \frac{1}{\theta^2\Delta t} \boldsymbol{v}_{n} - \frac{1-\theta}{\theta}\boldsymbol{a}_{n} \end{aligned}\end{split}\]

• option eval_nonlin_terms : "midpoint":

(16)\[\begin{split}\begin{aligned} \boldsymbol{u}_{n+\theta} &= \theta \boldsymbol{u}_{n+1} + (1-\theta) \boldsymbol{u}_{n} \\ \boldsymbol{v}_{n+\theta} &= \theta \boldsymbol{v}_{n+1} + (1-\theta) \boldsymbol{v}_{n} \\ \boldsymbol{a}_{n+\theta} &= \theta \boldsymbol{a}_{n+1} + (1-\theta) \boldsymbol{a}_{n} \end{aligned}\end{split}\]

\[\begin{aligned} \delta \mathcal{W}_{\mathrm{kin}}(\boldsymbol{a}_{n+\theta};\delta\boldsymbol{u}) + \delta \mathcal{W}_{\mathrm{int}}(\boldsymbol{u}_{n+\theta};\delta\boldsymbol{u}) - \delta \mathcal{W}_{\mathrm{ext}}(\boldsymbol{u}_{n+\theta};\delta\boldsymbol{u}) = 0, \quad \forall \; \delta\boldsymbol{u}\end{aligned}\]

• option eval_nonlin_terms : "trapezoidal":

\[\begin{split}\begin{aligned} &\theta\,\delta \mathcal{W}_{\mathrm{kin}}(\boldsymbol{a}_{n+1};\delta\boldsymbol{u}) + (1-\theta)\,\delta \mathcal{W}_{\mathrm{kin}}(\boldsymbol{a}_{n};\delta\boldsymbol{u}) + \\ & \theta\,\delta \mathcal{W}_{\mathrm{int}}(\boldsymbol{u}_{n+1};\delta\boldsymbol{u}) + (1-\theta)\,\delta \mathcal{W}_{\mathrm{int}}(\boldsymbol{u}_{n};\delta\boldsymbol{u}) - \\ & \theta\,\delta \mathcal{W}_{\mathrm{ext}}(\boldsymbol{u}_{n+1};\delta\boldsymbol{u}) - (1-\theta)\,\delta \mathcal{W}_{\mathrm{ext}}(\boldsymbol{u}_{n};\delta\boldsymbol{u}) = 0, \quad \forall \; \delta\boldsymbol{u} \end{aligned}\end{split}\]

Note the equivalence of "midpoint" and "trapezoidal" for all linear terms, e.g. \(\delta \mathcal{W}_{\mathrm{kin}}\), or for no or only linear dependence of \(\delta \mathcal{W}_{\mathrm{ext}}\) on the solution.

Note that, for incompressible mechanics, the pressure kinematic constraint is always evaluated at \(t_{n+1}\):

\[\begin{aligned} \delta \mathcal{W}_{\mathrm{pres}}(\boldsymbol{u}_{n+1};\delta p) = 0, \quad \forall \; \delta p, \end{aligned}\]

and the pressure in \(\delta \mathcal{W}_{\mathrm{int}}\) is set according to ([equation-solid-midpoint-genalpha]) or ([equation-solid-midpoint-ost]), respectively.

Spatial discretization and solution

– Discrete nonlinear system to solve in each time step \(n\) (displacement-based):

(17)\[\begin{aligned} \left.\boldsymbol{\mathsf{r}}_{u}(\boldsymbol{\mathsf{u}})\right|_{n+1} = \boldsymbol{\mathsf{0}} \end{aligned}\]

– Discrete linear system to solve in each Newton iteration \(k\) (displacement-based):

(18)\[\begin{aligned} \left. \boldsymbol{\mathsf{K}}_{uu} \right|_{n+1}^{k} \Delta\boldsymbol{\mathsf{u}}_{n+1}^{k+1}=-\left. \boldsymbol{\mathsf{r}}_{u} \right|_{n+1}^{k} \end{aligned}\]

– Discrete nonlinear system to solve in each time step \(n\) (incompressible):

(19)\[\begin{split}\begin{aligned} \boldsymbol{\mathsf{r}}_{n+1} = \begin{bmatrix} \boldsymbol{\mathsf{r}}_{u}(\boldsymbol{\mathsf{u}},\boldsymbol{\mathsf{p}}) \\ \boldsymbol{\mathsf{r}}_{p}(\boldsymbol{\mathsf{u}}) \end{bmatrix}_{n+1} = \boldsymbol{\mathsf{0}} \end{aligned}\end{split}\]

– Discrete linear system to solve in each Newton iteration \(k\) (incompressible):

(20)\[\begin{split}\begin{aligned} \begin{bmatrix} \boldsymbol{\mathsf{K}}_{uu} & \boldsymbol{\mathsf{K}}_{up} \\ \\ \boldsymbol{\mathsf{K}}_{pu} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}\end{bmatrix}_{n+1}^{k}\begin{bmatrix} \Delta\boldsymbol{\mathsf{u}} \\ \\ \Delta\boldsymbol{\mathsf{p}} \end{bmatrix}_{n+1}^{k+1}=-\begin{bmatrix} \boldsymbol{\mathsf{r}}_{u} \\ \\ \boldsymbol{\mathsf{r}}_{p} \end{bmatrix}_{n+1}^{k} \end{aligned}\end{split}\]

Fluid mechanics

Eulerian reference frame

– Example: Sec. 5.2 and demos/fluid

– Problem type: fluid

– Incompressible Navier-Stokes equations in Eulerian reference frame

Strong Form

– Primary variables: velocity \(\boldsymbol{v}\) and pressure \(p\)

(21)\[\begin{split}\begin{aligned} \rho\left(\frac{\partial\boldsymbol{v}}{\partial t} + (\nabla\boldsymbol{v})\,\boldsymbol{v}\right) &= \nabla \cdot \boldsymbol{\sigma}(\boldsymbol{v},p) + \hat{\boldsymbol{b}} &&\text{in} \; \mathit{\mathit{\Omega}}_t \times [0, T], \\ \nabla\cdot \boldsymbol{v} &= 0 &&\text{in} \; \mathit{\mathit{\Omega}}_t \times [0, T],\\ \boldsymbol{v} &= \hat{\boldsymbol{v}} &&\text{on} \; \mathit{\mathit{\Gamma}}_t^{\mathrm{D}} \times [0, T],\\ \boldsymbol{t} = \boldsymbol{\sigma}\boldsymbol{n} &= \hat{\boldsymbol{t}} &&\text{on} \; \mathit{\mathit{\Gamma}}_t^{\mathrm{N}} \times [0, T],\\ \boldsymbol{v}(\boldsymbol{x},0) &= \hat{\boldsymbol{v}}_{0}(\boldsymbol{x}) &&\text{in} \; \mathit{\mathit{\Omega}}_t, \end{aligned}\end{split}\]

with a Newtonian fluid constitutive law

\[\begin{aligned} \boldsymbol{\sigma} = -p \boldsymbol{I} + \eta \left(\nabla \boldsymbol{v} + (\nabla \boldsymbol{v})^{\mathrm{T}}\right) \end{aligned}\]

where \(\eta\) is the dynamic viscosity

Weak Form

Find \((\boldsymbol{v},p) \in \boldsymbol{\mathcal{V}}_{v}^{h,D} \times \mathcal{V}_{p}^{h,D}\) such that

(22)\[\begin{split}\begin{aligned} r_v(\boldsymbol{v},p;\delta\boldsymbol{v}) &:= \delta \mathcal{P}_{\mathrm{kin}}(\boldsymbol{v};\delta\boldsymbol{v}) + \delta \mathcal{P}_{\mathrm{int}}(\boldsymbol{v},p;\delta\boldsymbol{v}) - \delta \mathcal{P}_{\mathrm{ext}}(\boldsymbol{v};\delta\boldsymbol{v}) = 0, \quad \\ r_p(\boldsymbol{v};\delta p) &:= \delta \mathcal{P}_{\mathrm{pres}}(\boldsymbol{v};\delta p) = 0, \end{aligned}\end{split}\]

for all \((\delta\boldsymbol{v},\delta p) \in \boldsymbol{\mathcal{V}}_{v}^{h} \times \mathcal{V}_{p}^{h}\)

– Kinetic virtual power:

(23)\[\begin{aligned} \delta \mathcal{P}_{\mathrm{kin}}(\boldsymbol{v};\delta\boldsymbol{v}) = \int\limits_{\mathit{\Omega}_t} \rho\left(\frac{\partial\boldsymbol{v}}{\partial t} + (\nabla\boldsymbol{v})\,\boldsymbol{v}\right) \cdot \delta\boldsymbol{v} \,\mathrm{d}v \end{aligned}\]

– Internal virtual power:

(24)\[\begin{aligned} \delta \mathcal{P}_{\mathrm{int}}(\boldsymbol{v},p;\delta\boldsymbol{v}) = \int\limits_{\mathit{\Omega}_t} \boldsymbol{\sigma}(\boldsymbol{v},p) : \nabla \delta\boldsymbol{v} \,\mathrm{d}v \end{aligned}\]

– Pressure virtual power:

(25)\[\begin{aligned} \delta \mathcal{P}_{\mathrm{pres}}(\boldsymbol{v};\delta p) = \int\limits_{\mathit{\Omega}_t} (\nabla\cdot\boldsymbol{v})\,\delta p\,\mathrm{d}v \end{aligned}\]

– External virtual power:

• conservative Neumann load:

  (26)\[\begin{aligned} \delta \mathcal{P}_{\mathrm{ext}}(\delta\boldsymbol{v}) &= \int\limits_{\mathit{\Gamma}_t^{\mathrm{N}}} \hat{\boldsymbol{t}}(t) \cdot \delta\boldsymbol{v} \,\mathrm{d}a \end{aligned}\]

• pressure Neumann load:

  (27)\[\begin{aligned} \delta \mathcal{P}_{\mathrm{ext}}(\delta\boldsymbol{v}) &= -\int\limits_{\mathit{\Gamma}_t^{\mathrm{N}}} \hat{p}(t)\,\boldsymbol{n} \cdot \delta\boldsymbol{v} \,\mathrm{d}a \end{aligned}\]

• body force:

  (28)\[\begin{aligned} \delta \mathcal{P}_{\mathrm{ext}}(\delta\boldsymbol{v}) &= \int\limits_{\mathit{\Omega}_t} \hat{\boldsymbol{b}}(t) \cdot \delta\boldsymbol{v} \,\mathrm{d}V \end{aligned}\]

Stabilization

Streamline-upwind Petrov-Galerkin/pressure-stabilizing Petrov-Galerkin (SUPG/PSPG) methods are implemented, either using the full or a reduced scheme

– to the fluid FEM params, add the dict entry:

"stabilization" : {"scheme" : <SCHEME>, "vscale" : 1e3, "dscales" : [<d1>,<d2>,<d3>],
                   "symmetric" : False, "reduced_scheme" : False}

Full scheme according to [TO00]: "supg_pspg":

– Velocity residual operator in ([equation-fluid-weak-form]) is augmented by the following terms:

\[\begin{split}\begin{aligned} r_v \leftarrow r_v &+ \frac{1}{\rho}\int\limits_{\mathit{\Omega}_t} \tau_{\mathrm{SUPG}}\,(\nabla\delta\boldsymbol{v})\,\boldsymbol{v} \cdot \left(\rho\left(\frac{\partial \boldsymbol{v}}{\partial t} + (\nabla\boldsymbol{v})\,\boldsymbol{v}\right) - \nabla \cdot \boldsymbol{\sigma}(\boldsymbol{v},p)\right)\,\mathrm{d}v \\ & + \int\limits_{\mathit{\Omega}_t} \tau_{\mathrm{LSIC}}\,\rho\,(\nabla\cdot\delta\boldsymbol{v})(\nabla\cdot\boldsymbol{v})\,\mathrm{d}v \end{aligned}\end{split}\]

– Pressure residual operator in ([equation-fluid-weak-form]) is augmented by the following terms:

\[\begin{aligned} r_p \leftarrow r_p &+ \frac{1}{\rho}\int\limits_{\mathit{\Omega}_t} \tau_{\mathrm{PSPG}}\,(\nabla\delta p) \cdot \left(\rho\left(\frac{\partial \boldsymbol{v}}{\partial t} + (\nabla\boldsymbol{v})\,\boldsymbol{v}\right) - \nabla \cdot \boldsymbol{\sigma}(\boldsymbol{v},p)\right)\,\mathrm{d}v \end{aligned}\]

– Discrete nonlinear system to solve in each time step \(n\):

(29)\[\begin{split}\begin{aligned} \boldsymbol{\mathsf{r}}_{n+1} = \begin{bmatrix} \boldsymbol{\mathsf{r}}_{v}(\boldsymbol{\mathsf{v}},\boldsymbol{\mathsf{p}}) \\ \boldsymbol{\mathsf{r}}_{p}(\boldsymbol{\mathsf{p}},\boldsymbol{\mathsf{v}}) \end{bmatrix}_{n+1} = \boldsymbol{\mathsf{0}} \end{aligned}\end{split}\]

– Discrete linear system to solve in each Newton iteration \(k\):

(30)\[\begin{split}\begin{aligned} \begin{bmatrix} \boldsymbol{\mathsf{K}}_{vv} & \boldsymbol{\mathsf{K}}_{vp} \\ \\ \boldsymbol{\mathsf{K}}_{pv} & \boldsymbol{\mathsf{K}}_{pp} \end{bmatrix}_{n+1}^{k}\begin{bmatrix} \Delta\boldsymbol{\mathsf{v}} \\ \\ \Delta\boldsymbol{\mathsf{p}} \end{bmatrix}_{n+1}^{k+1}=-\begin{bmatrix} \boldsymbol{\mathsf{r}}_{v} \\ \\ \boldsymbol{\mathsf{r}}_{p} \end{bmatrix}_{n+1}^{k} \end{aligned}\end{split}\]

– Note that \(\boldsymbol{\mathsf{K}}_{pp}\) is zero for Taylor-Hood elements (without stabilization)

ALE reference frame

– Problem type: fluid_ale

– Incompressible Navier-Stokes equations in Arbitrary Lagrangian Eulerian (ALE) reference frame

– ALE domain problem deformation governed by linear-elastic, nonlinear hyperelastic solid, or a diffusion problem, displacement field \(\boldsymbol{d}\)

– Fluid mechanics formulated with respect to the reference frame, using ALE deformation gradient \(\widehat{\boldsymbol{F}}(\boldsymbol{d}) = \boldsymbol{I} + \nabla_0\boldsymbol{d}\) and its determinant, \(\widehat{J}(\boldsymbol{d})=\det \widehat{\boldsymbol{F}}(\boldsymbol{d})\)

ALE problem

– Primary variable: domain displacement \(\boldsymbol{d}\)

– Strong form:

(31)\[\begin{split}\begin{aligned} \nabla_{0} \cdot \boldsymbol{\sigma}^{\mathrm{G}}(\boldsymbol{d}) &= \boldsymbol{0} &&\text{in} \; \mathit{\mathit{\Omega}}_0, \\ \boldsymbol{d} &= \hat{\boldsymbol{d}} &&\text{on} \; \mathit{\mathit{\Gamma}}_0^{\mathrm{D}}, \end{aligned}\end{split}\]

– ALE material linelast:

\[\begin{aligned} \boldsymbol{\sigma}^{\mathrm{G}}(\boldsymbol{d}) = 2\mu \,\boldsymbol{\varepsilon} + \lambda \,\mathrm{tr}\boldsymbol{\varepsilon}\,\boldsymbol{I}, \qquad \text{with}\quad \boldsymbol{\varepsilon} = \frac{1}{2}\left(\nabla_0\boldsymbol{d} + (\nabla_0\boldsymbol{d})^{\mathrm{T}}\right) \end{aligned}\]

– ALE material diffusion:

\[\begin{aligned} \boldsymbol{\sigma}^{\mathrm{G}}(\boldsymbol{d}) = D \,\nabla_0\boldsymbol{d} \end{aligned}\]

– ALE material neohooke (fully nonlinear model):

\[\begin{aligned} \boldsymbol{\sigma}^{\mathrm{G}}(\boldsymbol{d}) = \frac{\partial \mathit{\Psi}}{\partial \widehat{\boldsymbol{F}}}, \qquad \text{with}\quad \mathit{\Psi} = \frac{\mu}{2}\left(\mathrm{tr}(\widehat{\boldsymbol{F}}^{\mathrm{T}}\widehat{\boldsymbol{F}}) - 3\right) + \frac{\mu}{2\beta} \left(\widehat{J}^{-2\beta} - 1\right) \end{aligned}\]

– weak form:

(32)\[\begin{aligned} r_{d}(\boldsymbol{d};\delta\boldsymbol{d}) := \int\limits_{\mathit{\Omega}_0}\boldsymbol{\sigma}^{\mathrm{G}}(\boldsymbol{d}) : \nabla_{0}\delta\boldsymbol{d}\,\mathrm{d}V = 0, \quad \forall \; \delta\boldsymbol{d} \end{aligned}\]

Strong form (ALE)

– Primary variables: velocity \(\boldsymbol{v}\), pressure \(p\), and domain displacement \(\boldsymbol{d}\)

(33)\[\begin{split}\begin{aligned} \widehat{J}\rho\left(\left.\frac{\partial\boldsymbol{v}}{\partial t}\right|_{\boldsymbol{x}_{0}} + (\nabla_0\boldsymbol{v}\,\widehat{\boldsymbol{F}}^{-1})\,(\boldsymbol{v}-\widehat{\boldsymbol{w}})\right) &= \nabla_{0} \cdot\left(\widehat{J}\boldsymbol{\sigma} \widehat{\boldsymbol{F}}^{-\mathrm{T}}\right) + \widehat{J}\hat{\boldsymbol{b}} &&\text{in} \; \mathit{\mathit{\Omega}}_0 \times [0, T],\\ \nabla_{0}\cdot\left(\widehat{J}\widehat{\boldsymbol{F}}^{-1}\boldsymbol{v}\right) &= 0 &&\text{in} \; \mathit{\mathit{\Omega}}_0 \times [0, T],\\ \boldsymbol{v} &= \hat{\boldsymbol{v}} &&\text{on} \; \mathit{\mathit{\Gamma}}_0^{\mathrm{D}} \times [0, T], \\ \boldsymbol{t} = \boldsymbol{\sigma}\boldsymbol{n} &= \hat{\boldsymbol{t}} &&\text{on} \; \mathit{\mathit{\Gamma}}_0^{\mathrm{N}} \times [0, T], \\ \boldsymbol{v}(\boldsymbol{x},0) &= \hat{\boldsymbol{v}}_{0}(\boldsymbol{x}) &&\text{in} \; \mathit{\mathit{\Omega}}_0, \end{aligned}\end{split}\]

with a Newtonian fluid constitutive law

\[\begin{aligned} \boldsymbol{\sigma} = -p \boldsymbol{I} + \mu \left(\nabla_0 \boldsymbol{v}\,\widehat{\boldsymbol{F}}^{-1} + \widehat{\boldsymbol{F}}^{-\mathrm{T}}(\nabla_0 \boldsymbol{v})^{\mathrm{T}}\right) \end{aligned}\]

Weak form (ALE)

Find \((\boldsymbol{v},p,\boldsymbol{d}) \in \boldsymbol{\mathcal{V}}_{v}^{h,D} \times \mathcal{V}_{p}^{h,D} \times \boldsymbol{\mathcal{V}}_{d}^{h,D}\) such that

(34)\[\begin{split}\begin{aligned} r_v(\boldsymbol{v},p,\boldsymbol{d};\delta\boldsymbol{v}) &:= \delta \mathcal{P}_{\mathrm{kin}}(\boldsymbol{v},\boldsymbol{d};\delta\boldsymbol{v}) + \delta \mathcal{P}_{\mathrm{int}}(\boldsymbol{v},p,\boldsymbol{d};\delta\boldsymbol{v}) - \delta \mathcal{P}_{\mathrm{ext}}(\boldsymbol{v},\boldsymbol{d};\delta\boldsymbol{v}) = 0, \\ r_p(\boldsymbol{v},\boldsymbol{d};\delta p) &:= \delta \mathcal{P}_{\mathrm{pres}}(\boldsymbol{v},\boldsymbol{d};\delta p) = 0, \\ r_{d}(\boldsymbol{d};\delta\boldsymbol{d}) &= 0, \end{aligned}\end{split}\]

for all \((\delta\boldsymbol{v},\delta p,\delta\boldsymbol{d}) \in \boldsymbol{\mathcal{V}}_{v}^{h} \times \mathcal{V}_{p}^{h} \times \boldsymbol{\mathcal{V}}_{d}^{h}\)

– Kinetic virtual power:

\[\begin{aligned} \delta \mathcal{P}_{\mathrm{kin}}(\boldsymbol{v},\boldsymbol{d};\delta\boldsymbol{v}) = \int\limits_{\mathit{\Omega}_0} \widehat{J} \rho\left(\left.\frac{\partial\boldsymbol{v}}{\partial t}\right|_{\boldsymbol{x}_{0}} + (\nabla_{0}\boldsymbol{v}\,\widehat{\boldsymbol{F}}^{-1})\,(\boldsymbol{v}-\widehat{\boldsymbol{w}})\right) \cdot \delta\boldsymbol{v} \,\mathrm{d}V \end{aligned}\]

– Internal virtual power:

\[\begin{aligned} \delta \mathcal{P}_{\mathrm{int}}(\boldsymbol{v},p,\boldsymbol{d};\delta\boldsymbol{v}) = \int\limits_{\mathit{\Omega}_0} \widehat{J}\boldsymbol{\sigma}\widehat{\boldsymbol{F}}^{-\mathrm{T}} : \nabla_{0} \delta\boldsymbol{v} \,\mathrm{d}V \end{aligned}\]

– Pressure virtual power:

\[\begin{aligned} \delta \mathcal{P}_{\mathrm{pres}}(\boldsymbol{v},\boldsymbol{d};\delta p) = \int\limits_{\mathit{\Omega}_0} \nabla_{0}\cdot\left(\widehat{J}\widehat{\boldsymbol{F}}^{-1}\boldsymbol{v}\right) \delta p\,\mathrm{d}V \end{aligned}\]

– External virtual power:

• conservative Neumann load:

  \[\begin{aligned} \delta \mathcal{P}_{\mathrm{ext}}(\delta\boldsymbol{v}) &= \int\limits_{\mathit{\Gamma}_0^{\mathrm{N}}} \hat{\boldsymbol{t}}(t) \cdot \delta\boldsymbol{v} \,\mathrm{d}A \end{aligned}\]

• pressure Neumann load:

  \[\begin{aligned} \delta \mathcal{P}_{\mathrm{ext}}(\boldsymbol{d};\delta\boldsymbol{v}) &= -\int\limits_{\mathit{\Gamma}_0^{\mathrm{N}}} \hat{p}(t)\,\widehat{J}\widehat{\boldsymbol{F}}^{-\mathrm{T}}\boldsymbol{n}_{0} \cdot \delta\boldsymbol{v} \,\mathrm{d}A \end{aligned}\]

• body force:

  \[\begin{aligned} \delta \mathcal{P}_{\mathrm{ext}}(\boldsymbol{d};\delta\boldsymbol{v}) &= \int\limits_{\mathit{\Omega}_0} \widehat{J}\,\hat{\boldsymbol{b}}(t) \cdot \delta\boldsymbol{v} \,\mathrm{d}V \end{aligned}\]

Stabilization (ALE)

"supg_pspg":

– Velocity residual operator in ([equation-fluid-ale-weak-form]) is augmented by the following terms:

\[\begin{split}\begin{aligned} r_v \leftarrow r_v &+ \frac{1}{\rho}\int\limits_{\mathit{\Omega}_0}\widehat{J}\, \tau_{\mathrm{SUPG}}\,(\nabla_0\delta\boldsymbol{v}\,\widehat{\boldsymbol{F}}^{-1})\,\boldsymbol{v}\;\cdot \\ & \qquad\quad \cdot\left(\widehat{J}\rho\left(\left.\frac{\partial \boldsymbol{v}}{\partial t}\right|_{\boldsymbol{x}_{0}} + (\nabla_0\boldsymbol{v}\,\widehat{\boldsymbol{F}}^{-1})\,(\boldsymbol{v}-\widehat{\boldsymbol{w}})\right) - \nabla_{0}\cdot\left(\widehat{J}\boldsymbol{\sigma}\widehat{\boldsymbol{F}}^{-\mathrm{T}}\right)\right)\,\mathrm{d}V \\ & + \int\limits_{\mathit{\Omega}_0}\widehat{J}\, \tau_{\mathrm{LSIC}}\,\rho\,(\nabla_{0}\delta\boldsymbol{v} : \widehat{\boldsymbol{F}}^{-\mathrm{T}})(\nabla_{0}\boldsymbol{v} : \widehat{\boldsymbol{F}}^{-\mathrm{T}})\,\mathrm{d}V \end{aligned}\end{split}\]

– Pressure residual operator in ([equation-fluid-ale-weak-form]) is augmented by the following terms:

\[\begin{split}\begin{aligned} r_p \leftarrow r_p &+ \frac{1}{\rho}\int\limits_{\mathit{\Omega}_0}\widehat{J}\, \tau_{\mathrm{PSPG}}\,(\widehat{\boldsymbol{F}}^{-\mathrm{T}}\nabla_{0}\delta p) \;\cdot \\ & \qquad\quad \cdot \left(\widehat{J}\rho\left(\left.\frac{\partial \boldsymbol{v}}{\partial t}\right|_{\boldsymbol{x}_{0}} + (\nabla_0\boldsymbol{v}\,\widehat{\boldsymbol{F}}^{-1})\,(\boldsymbol{v}-\widehat{\boldsymbol{w}})\right) - \nabla_{0}\cdot\left(\widehat{J}\boldsymbol{\sigma}\widehat{\boldsymbol{F}}^{-\mathrm{T}}\right)\right)\,\mathrm{d}V \end{aligned}\end{split}\]

– Discrete nonlinear system to solve in each time step \(n\):

(35)\[\begin{split}\begin{aligned} \boldsymbol{\mathsf{r}}_{n+1} = \begin{bmatrix} \boldsymbol{\mathsf{r}}_{v}(\boldsymbol{\mathsf{v}},\boldsymbol{\mathsf{p}},\boldsymbol{\mathsf{d}}) \\ \boldsymbol{\mathsf{r}}_{p}(\boldsymbol{\mathsf{p}},\boldsymbol{\mathsf{v}},\boldsymbol{\mathsf{d}}) \\ \boldsymbol{\mathsf{r}}_{d}(\boldsymbol{\mathsf{d}}) \end{bmatrix}_{n+1} = \boldsymbol{\mathsf{0}} \end{aligned}\end{split}\]

– Discrete linear system to solve in each Newton iteration \(k\):

(36)\[\begin{split}\begin{aligned} \begin{bmatrix} \boldsymbol{\mathsf{K}}_{vv} & \boldsymbol{\mathsf{K}}_{vp} & \boldsymbol{\mathsf{K}}_{vd} \\ \\ \boldsymbol{\mathsf{K}}_{pv} & \boldsymbol{\mathsf{K}}_{pp} & \boldsymbol{\mathsf{K}}_{pd} \\ \\ \boldsymbol{\mathsf{K}}_{dv} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{dd} \end{bmatrix}_{n+1}^{k}\begin{bmatrix} \Delta\boldsymbol{\mathsf{v}} \\ \\ \Delta\boldsymbol{\mathsf{p}} \\ \\ \Delta\boldsymbol{\mathsf{d}} \end{bmatrix}_{n+1}^{k+1}=-\begin{bmatrix} \boldsymbol{\mathsf{r}}_{v} \\ \\ \boldsymbol{\mathsf{r}}_{p} \\ \\ \boldsymbol{\mathsf{r}}_{d} \end{bmatrix}_{n+1}^{k} \end{aligned}\end{split}\]

– note that \(\boldsymbol{\mathsf{K}}_{pp}\) is zero for Taylor-Hood elements (without stabilization)

0D flow: Lumped parameter models

– Example: Sec. 5.3 and demos/flow0d

– Problem type: flow0d

– 0D model concentrated elements are resistances (\(R\)), impedances (\(Z\), technically are resistances), compliances (\(C\)), inertances (\(L\) or \(I\)), and elastances (\(E\))

– 0D variables are pressures (\(p\)), fluxes (\(q\) or \(Q\)), or volumes (\(V\))

2-element Windkessel

– Model type : 2elwindkessel

4-element Windkessel (inertance parallel to impedance)

– Model type : 4elwindkesselLpZ

4-element Windkessel (inertance serial to impedance)

– Model type : 4elwindkesselLsZ

In-outflow CRL link

– Model type : CRLinoutlink

Systemic and pulmonary circulation

– Model type : syspul

– Allows to link in a coronary flow model

(37)\[\begin{split}\begin{aligned} &\text{left heart and systemic circulation} && \nonumber\\ &-Q_{\mathrm{at}}^{\ell} = \sum\limits_{i=1}^{n_{\mathrm{ven}}^{\mathrm{pul}}}q_{\mathrm{ven},i}^{\mathrm{pul}} - q_{\mathrm{v,in}}^{\ell} && \text{left atrium flow balance}\\ &q_{\mathrm{v,in}}^{\ell} = q_{\mathrm{mv}}(p_{\mathrm{at}}^{\ell}-p_{\mathrm{v}}^{\ell}) && \text{mitral valve momentum}\\ &-Q_{\mathrm{v}}^{\ell} = q_{\mathrm{v,in}}^{\ell} - q_{\mathrm{v,out}}^{\ell} && \text{left ventricle flow balance}\\ &q_{\mathrm{v,out}}^{\ell} = q_{\mathrm{av}}(p_{\mathrm{v}}^{\ell}-p_{\mathrm{ar}}^{\mathrm{sys}}) && \text{aortic valve momentum}\\ &-Q_{\mathrm{aort}}^{\mathrm{sys}} = q_{\mathrm{v,out}}^{\ell} - q_{\mathrm{ar,p}}^{\mathrm{sys}} - \mathbb{I}^{\mathrm{cor}}\sum\limits_{i\in\{\ell,r\}}q_{\mathrm{cor,p,in}}^{\mathrm{sys},i} && \text{aortic root flow balance}\\ &I_{\mathrm{ar}}^{\mathrm{sys}} \frac{\mathrm{d}q_{\mathrm{ar,p}}^{\mathrm{sys}}}{\mathrm{d}t} + Z_{\mathrm{ar}}^{\mathrm{sys}}\,q_{\mathrm{ar,p}}^{\mathrm{sys}}=p_{\mathrm{ar}}^{\mathrm{sys}}-p_{\mathrm{ar,d}}^{\mathrm{sys}} && \text{aortic root inertia}\nonumber\\ &C_{\mathrm{ar}}^{\mathrm{sys}} \frac{\mathrm{d}p_{\mathrm{ar,d}}^{\mathrm{sys}}}{\mathrm{d}t} = q_{\mathrm{ar,p}}^{\mathrm{sys}} - q_{\mathrm{ar}}^{\mathrm{sys}} && \text{systemic arterial flow balance}\\ &L_{\mathrm{ar}}^{\mathrm{sys}} \frac{\mathrm{d}q_{\mathrm{ar}}^{\mathrm{sys}}}{\mathrm{d}t} + R_{\mathrm{ar}}^{\mathrm{sys}}\,q_{\mathrm{ar}}^{\mathrm{sys}}=p_{\mathrm{ar,d}}^{\mathrm{sys}}-p_{\mathrm{ven}}^{\mathrm{sys}} && \text{systemic arterial momentum}\\ &C_{\mathrm{ven}}^{\mathrm{sys}} \frac{\mathrm{d}p_{\mathrm{ven}}^{\mathrm{sys}}}{\mathrm{d}t} = q_{\mathrm{ar}}^{\mathrm{sys}}-\sum\limits_{i=1}^{n_{\mathrm{ven}}^{\mathrm{sys}}}q_{\mathrm{ven},i}^{\mathrm{sys}}\ && \text{systemic venous flow balance}\\ &L_{\mathrm{ven},i}^{\mathrm{sys}} \frac{\mathrm{d}q_{\mathrm{ven},i}^{\mathrm{sys}}}{\mathrm{d}t} + R_{\mathrm{ven},i}^{\mathrm{sys}}\, q_{\mathrm{ven},i}^{\mathrm{sys}} = p_{\mathrm{ven}}^{\mathrm{sys}} - p_{\mathrm{at},i}^{r} && \text{systemic venous momentum}\nonumber\\ &\qquad\qquad i \in \{1,...,n_{\mathrm{ven}}^{\mathrm{sys}}\} && \end{aligned}\end{split}\]

(38)\[\begin{split}\begin{aligned} &\text{right heart and pulmonary circulation} && \\ &-Q_{\mathrm{at}}^{r} = \sum\limits_{i=1}^{n_{\mathrm{ven}}^{\mathrm{sys}}}q_{\mathrm{ven},i}^{\mathrm{sys}} + \mathbb{I}^{\mathrm{cor}} q_{\mathrm{cor,d,out}}^{\mathrm{sys}} - q_{\mathrm{v,in}}^{r} && \text{right atrium flow balance}\\ &q_{\mathrm{v,in}}^{r} = q_{\mathrm{tv}}(p_{\mathrm{at}}^{r}-p_{\mathrm{v}}^{r}) && \text{tricuspid valve momentum}\\ &-Q_{\mathrm{v}}^{r} = q_{\mathrm{v,in}}^{r} - q_{\mathrm{v,out}}^{r} && \text{right ventricle flow balance}\\ &q_{\mathrm{v,out}}^{r} = q_{\mathrm{pv}}(p_{\mathrm{v}}^{r}-p_{\mathrm{ar}}^{\mathrm{pul}}) && \text{pulmonary valve momentum}\\ &C_{\mathrm{ar}}^{\mathrm{pul}} \frac{\mathrm{d}p_{\mathrm{ar}}^{\mathrm{pul}}}{\mathrm{d}t} = q_{\mathrm{v,out}}^{r} - q_{\mathrm{ar}}^{\mathrm{pul}} && \text{pulmonary arterial flow balance}\\ &L_{\mathrm{ar}}^{\mathrm{pul}} \frac{\mathrm{d}q_{\mathrm{ar}}^{\mathrm{pul}}}{\mathrm{d}t} + R_{\mathrm{ar}}^{\mathrm{pul}}\,q_{\mathrm{ar}}^{\mathrm{pul}}=p_{\mathrm{ar}}^{\mathrm{pul}} -p_{\mathrm{ven}}^{\mathrm{pul}} && \text{pulmonary arterial momentum}\\ &C_{\mathrm{ven}}^{\mathrm{pul}} \frac{\mathrm{d}p_{\mathrm{ven}}^{\mathrm{pul}}}{\mathrm{d}t} = q_{\mathrm{ar}}^{\mathrm{pul}} - \sum\limits_{i=1}^{n_{\mathrm{ven}}^{\mathrm{pul}}}q_{\mathrm{ven},i}^{\mathrm{pul}} && \text{pulmonary venous flow balance}\nonumber\\ &L_{\mathrm{ven},i}^{\mathrm{pul}} \frac{\mathrm{d}q_{\mathrm{ven},i}^{\mathrm{pul}}}{\mathrm{d}t} + R_{\mathrm{ven},i}^{\mathrm{pul}}\, q_{\mathrm{ven},i}^{\mathrm{pul}}=p_{\mathrm{ven}}^{\mathrm{pul}}-p_{\mathrm{at},i}^{\ell} && \text{pulmonary venous momentum}\\ &\qquad\qquad i \in \{1,...,n_{\mathrm{ven}}^{\mathrm{pul}}\} && \end{aligned}\end{split}\]

with:

\[\begin{aligned} Q_{\mathrm{at}}^{\ell} := -\frac{\mathrm{d}V_{\mathrm{at}}^{\ell}}{\mathrm{d}t}, \quad Q_{\mathrm{v}}^{\ell} := -\frac{\mathrm{d}V_{\mathrm{v}}^{\ell}}{\mathrm{d}t}, \quad Q_{\mathrm{at}}^{r} := -\frac{\mathrm{d}V_{\mathrm{at}}^{r}}{\mathrm{d}t}, \quad Q_{\mathrm{v}}^{r} := -\frac{\mathrm{d}V_{\mathrm{v}}^{r}}{\mathrm{d}t}, \quad Q_{\mathrm{aort}}^{\mathrm{sys}} := -\frac{\mathrm{d}V_{\mathrm{aort}}^{\mathrm{sys}}}{\mathrm{d}t} \end{aligned}\]

and:

\[\begin{split}\begin{aligned} \mathbb{I}^{\mathrm{cor}} = \begin{cases} 1, & \text{if \, CORONARYMODEL}, \\ 0, & \text{else} \end{cases} \end{aligned}\end{split}\]

The volume \(V\) of the heart chambers (0D) is modeled by the volume-pressure relationship

\[\begin{aligned} V(t) = \frac{p}{E(t)} + V_{\mathrm{u}}, \end{aligned}\]

with the unstressed volume \(V_{\mathrm{u}}\) and the time-varying elastance

(39)\[\begin{aligned} E(t)=\left(E_{\mathrm{max}}-E_{\mathrm{min}}\right)\cdot \hat{y}(t)+E_{\mathrm{min}}, \end{aligned}\]

where \(E_{\mathrm{max}}\) and \(E_{\mathrm{min}}\) denote the maximum and minimum elastance, respectively. The normalized activation function \(\hat{y}(t)\) is input by the user.

Flow-pressure relations for the four valves, eq. ([eq:mv_flow]), ([eq:av_flow]), ([eq:tv_flow]), ([eq:pv_flow]), are functions of the pressure difference \(p-p_{\mathrm{open}}\) across the valve. The following valve models can be defined:

Valve model pwlin_pres:

\[\begin{split}\begin{aligned} q(p-p_{\mathrm{open}}) = \frac{p-p_{\mathrm{open}}}{\tilde{R}}, \quad \text{with}\; \tilde{R} = \begin{cases} R_{\max}, & p < p_{\mathrm{open}} \\ R_{\min}, & p \geq p_{\mathrm{open}} \end{cases} \end{aligned}\end{split}\]

Remark: Non-smooth flow-pressure relationship

Valve model pwlin_time:

\[\begin{split}\begin{aligned} q(p-p_{\mathrm{open}}) = \frac{p-p_{\mathrm{open}}}{\tilde{R}},\quad \text{with}\; \tilde{R} = \begin{cases} \begin{cases} R_{\max}, & t < t_{\mathrm{open}} \;\text{and}\; t \geq t_{\mathrm{close}} \\ R_{\min}, & t \geq t_{\mathrm{open}} \;\text{or}\; t < t_{\mathrm{close}} \end{cases}, & t_{\mathrm{open}} > t_{\mathrm{close}} \\ \begin{cases} R_{\max}, & t < t_{\mathrm{open}} \;\text{or}\; t \geq t_{\mathrm{close}} \\ R_{\min}, & t \geq t_{\mathrm{open}} \;\text{and}\; t < t_{\mathrm{close}} \end{cases}, & \text{else} \end{cases} \end{aligned}\end{split}\]

Remark: Non-smooth flow-pressure relationship with resistance only dependent on timings, not the pressure difference!

Valve model smooth_pres_resistance:

\[\begin{aligned} q(p-p_{\mathrm{open}}) = \frac{p-p_{\mathrm{open}}}{\tilde{R}},\quad \text{with}\;\tilde{R} = 0.5\left(R_{\max}-R_{\min}\right)\left(\tanh\frac{p-p_{\mathrm{open}}}{\epsilon}+1\right) + R_{\min} \end{aligned}\]

Remark: Smooth but potentially non-convex flow-pressure relationship!

Valve model smooth_pres_momentum:

\[\begin{split}\begin{aligned} q(p-p_{\mathrm{open}}) = \begin{cases}\frac{p-p_{\mathrm{open}}}{R_{\max}}, & p < p_{\mathrm{open}}-0.5\epsilon \\ h_{00}p_{0} + h_{10}m_{0}\epsilon + h_{01}p_{1} + h_{11}m_{1}\epsilon, & p \geq p_{\mathrm{open}}-0.5\epsilon \;\text{and}\; p < p_{\mathrm{open}}+0.5\epsilon \\ \frac{p-p_{\mathrm{open}}}{R_{\min}}, & p \geq p_{\mathrm{open}}+0.5\epsilon \end{cases} \end{aligned}\end{split}\]

with

\[\begin{aligned} p_{0}=\frac{-0.5\epsilon}{R_{\max}}, \qquad m_{0}=\frac{1}{R_{\max}}, \qquad && p_{1}=\frac{0.5\epsilon}{R_{\min}}, \qquad m_{1}=\frac{1}{R_{\min}} \end{aligned}\]

and

\[\begin{split}\begin{aligned} h_{00}=2s^3 - 3s^2 + 1, &\qquad h_{01}=-2s^3 + 3s^2, \\ h_{10}=s^3 - 2s^2 + s, &\qquad h_{11}=s^3 - s^2 \end{aligned}\end{split}\]

with

\[\begin{aligned} s=\frac{p-p_{\mathrm{open}}+0.5\epsilon}{\epsilon} \end{aligned}\]

Remarks:

– Collapses to valve model pwlin_pres for \(\epsilon=0\)

– Smooth and convex flow-pressure relationship

Valve model pw_pres_regurg:

\[\begin{split}\begin{aligned} q(p-p_{\mathrm{open}}) = \begin{cases} c A_{\mathrm{o}} \sqrt{p-p_{\mathrm{open}}}, & p < p_{\mathrm{open}} \\ \frac{p-p_{\mathrm{open}}}{R_{\min}}, & p \geq p_{\mathrm{open}} \end{cases} \end{aligned}\end{split}\]

Remark: Model to allow a regurgitant valve in the closed state, degree of regurgitation can be varied by specifying the valve regurgitant area \(A_{\mathrm{o}}\)

– Coronary circulation model: ZCRp_CRd_lr

\[\begin{split}\begin{aligned} &C_{\mathrm{cor,p}}^{\mathrm{sys},\ell} \left(\frac{\mathrm{d}p_{\mathrm{ar}}^{\mathrm{sys},\ell}}{\mathrm{d}t}-Z_{\mathrm{cor,p}}^{\mathrm{sys},\ell}\frac{\mathrm{d}q_{\mathrm{cor,p,in}}^{\mathrm{sys},\ell}}{\mathrm{d}t}\right) = q_{\mathrm{cor,p,in}}^{\mathrm{sys},\ell} - q_{\mathrm{cor,p}}^{\mathrm{sys},\ell} && \text{left coronary proximal flow balance}\\ &R_{\mathrm{cor,p}}^{\mathrm{sys},\ell}\,q_{\mathrm{cor,p}}^{\mathrm{sys},\ell}=p_{\mathrm{ar}}^{\mathrm{sys}}-p_{\mathrm{cor,d}}^{\mathrm{sys},\ell} - Z_{\mathrm{cor,p}}^{\mathrm{sys},\ell}\,q_{\mathrm{cor,p,in}}^{\mathrm{sys},\ell} && \text{left coronary proximal momentum}\\ &C_{\mathrm{cor,d}}^{\mathrm{sys},\ell} \frac{\mathrm{d}(p_{\mathrm{cor,d}}^{\mathrm{sys},\ell}-p_{\mathrm{v}}^{\ell})}{\mathrm{d}t} = q_{\mathrm{cor,p}}^{\mathrm{sys},\ell} - q_{\mathrm{cor,d}}^{\mathrm{sys},\ell} && \text{left coronary distal flow balance}\\ &R_{\mathrm{cor,d}}^{\mathrm{sys},\ell}\,q_{\mathrm{cor,d}}^{\mathrm{sys},\ell}=p_{\mathrm{cor,d}}^{\mathrm{sys},\ell}-p_{\mathrm{at}}^{r} && \text{left coronary distal momentum}\\ &C_{\mathrm{cor,p}}^{\mathrm{sys},r} \left(\frac{\mathrm{d}p_{\mathrm{ar}}^{\mathrm{sys},r}}{\mathrm{d}t}-Z_{\mathrm{cor,p}}^{\mathrm{sys},r}\frac{\mathrm{d}q_{\mathrm{cor,p,in}}^{\mathrm{sys},r}}{\mathrm{d}t}\right) = q_{\mathrm{cor,p,in}}^{\mathrm{sys},r} - q_{\mathrm{cor,p}}^{\mathrm{sys},r} && \text{right coronary proximal flow balance}\\ &R_{\mathrm{cor,p}}^{\mathrm{sys},r}\,q_{\mathrm{cor,p}}^{\mathrm{sys},r}=p_{\mathrm{ar}}^{\mathrm{sys}}-p_{\mathrm{cor,d}}^{\mathrm{sys},r} - Z_{\mathrm{cor,p}}^{\mathrm{sys},r}\,q_{\mathrm{cor,p,in}}^{\mathrm{sys},r} && \text{right coronary proximal momentum}\\ &C_{\mathrm{cor,d}}^{\mathrm{sys},r} \frac{\mathrm{d}(p_{\mathrm{cor,d}}^{\mathrm{sys},r}-p_{\mathrm{v}}^{\ell})}{\mathrm{d}t} = q_{\mathrm{cor,p}}^{\mathrm{sys},r} - q_{\mathrm{cor,d}}^{\mathrm{sys},r} && \text{right coronary distal flow balance}\nonumber\\ &R_{\mathrm{cor,d}}^{\mathrm{sys},r}\,q_{\mathrm{cor,d}}^{\mathrm{sys},r}=p_{\mathrm{cor,d}}^{\mathrm{sys},r}-p_{\mathrm{at}}^{r} && \text{right coronary distal momentum}\\ &0=q_{\mathrm{cor,d}}^{\mathrm{sys},\ell}+q_{\mathrm{cor,d}}^{\mathrm{sys},r}-q_{\mathrm{cor,d,out}}^{\mathrm{sys}} && \text{distal coronary junction flow balance} \end{aligned}\end{split}\]

– Coronary circulation model: ZCRp_CRd

\[\begin{split}\begin{aligned} &C_{\mathrm{cor,p}}^{\mathrm{sys}} \left(\frac{\mathrm{d}p_{\mathrm{ar}}^{\mathrm{sys}}}{\mathrm{d}t}-Z_{\mathrm{cor,p}}^{\mathrm{sys}}\frac{\mathrm{d}q_{\mathrm{cor,p,in}}^{\mathrm{sys}}}{\mathrm{d}t}\right) = q_{\mathrm{cor,p,in}}^{\mathrm{sys}} - q_{\mathrm{cor,p}}^{\mathrm{sys}} && \text{coronary proximal flow balance}\\ &R_{\mathrm{cor,p}}^{\mathrm{sys}}\,q_{\mathrm{cor,p}}^{\mathrm{sys}}=p_{\mathrm{ar}}^{\mathrm{sys}}-p_{\mathrm{cor,d}}^{\mathrm{sys}} - Z_{\mathrm{cor,p}}^{\mathrm{sys}}\,q_{\mathrm{cor,p,in}}^{\mathrm{sys}} && \text{coronary proximal momentum}\nonumber\\ &C_{\mathrm{cor,d}}^{\mathrm{sys}} \frac{\mathrm{d}(p_{\mathrm{cor,d}}^{\mathrm{sys}}-p_{\mathrm{v}}^{\ell})}{\mathrm{d}t} = q_{\mathrm{cor,p}}^{\mathrm{sys}} - q_{\mathrm{cor,d}}^{\mathrm{sys}} && \text{coronary distal flow balance}\\ &R_{\mathrm{cor,d}}^{\mathrm{sys}}\,q_{\mathrm{cor,d}}^{\mathrm{sys}}=p_{\mathrm{cor,d}}^{\mathrm{sys}}-p_{\mathrm{at}}^{r} && \text{coronary distal momentum} \end{aligned}\end{split}\]

Systemic and pulmonary circulation, including capillary flow

– Model type : syspulcap, cf. [Hir19], p. 51ff.

(40)\[\begin{split}\begin{aligned} &-Q_{\mathrm{at}}^{\ell} = q_{\mathrm{ven}}^{\mathrm{pul}} - q_{\mathrm{v,in}}^{\ell} && \text{left atrium flow balance}\nonumber\\ &\tilde{R}_{\mathrm{v,in}}^{\ell}\,q_{\mathrm{v,in}}^{\ell} = p_{\mathrm{at}}^{\ell}-p_{\mathrm{v}}^{\ell} && \text{mitral valve momentum}\nonumber\\ &-Q_{\mathrm{v}}^{\ell} = q_{\mathrm{v,in}}^{\ell} - q_{\mathrm{v,out}}^{\ell} && \text{left ventricle flow balance}\nonumber\\ &\tilde{R}_{\mathrm{v,out}}^{\ell}\,q_{\mathrm{v,out}}^{\ell} = p_{\mathrm{v}}^{\ell}-p_{\mathrm{ar}}^{\mathrm{sys}} && \text{aortic valve momentum}\nonumber\\ &0 = q_{\mathrm{v,out}}^{\ell} - q_{\mathrm{ar,p}}^{\mathrm{sys}} && \text{aortic root flow balance}\nonumber\\ &I_{\mathrm{ar}}^{\mathrm{sys}} \frac{\mathrm{d}q_{\mathrm{ar,p}}^{\mathrm{sys}}}{\mathrm{d}t} + Z_{\mathrm{ar}}^{\mathrm{sys}}\,q_{\mathrm{ar,p}}^{\mathrm{sys}}=p_{\mathrm{ar}}^{\mathrm{sys}}-p_{\mathrm{ar,d}}^{\mathrm{sys}} && \text{aortic root inertia}\nonumber\\ &C_{\mathrm{ar}}^{\mathrm{sys}} \frac{\mathrm{d}p_{\mathrm{ar,d}}^{\mathrm{sys}}}{\mathrm{d}t} = q_{\mathrm{ar,p}}^{\mathrm{sys}} - q_{\mathrm{ar}}^{\mathrm{sys}} && \text{systemic arterial flow balance}\nonumber\\ &L_{\mathrm{ar}}^{\mathrm{sys}}\frac{\mathrm{d}q_{\mathrm{ar}}^{\mathrm{sys}}}{\mathrm{d}t} + R_{\mathrm{ar}}^{\mathrm{sys}}\,q_{\mathrm{ar}}^{\mathrm{sys}}=p_{\mathrm{ar,d}}^{\mathrm{sys}} -p_{\mathrm{ar,peri}}^{\mathrm{sys}} && \text{systemic arterial momentum}\nonumber\\ &\left(\sum_{j\in\{\mathrm{spl,espl,\atop msc,cer,cor}\}}\!\!\!\!\!\!\!\!\!C_{\mathrm{ar},j}^{\mathrm{sys}}\right) \frac{\mathrm{d}p_{\mathrm{ar,peri}}^{\mathrm{sys}}}{\mathrm{d}t} = q_{\mathrm{ar}}^{\mathrm{sys}}-\!\!\!\!\!\sum_{j\in\{\mathrm{spl,espl,\atop msc,cer,cor}\}}\!\!\!\!\!\!\!\!\!q_{\mathrm{ar},j}^{\mathrm{sys}} && \text{systemic capillary arterial flow balance}\nonumber\\ &R_{\mathrm{ar},i}^{\mathrm{sys}}\,q_{\mathrm{ar},i}^{\mathrm{sys}} = p_{\mathrm{ar,peri}}^{\mathrm{sys}} - p_{\mathrm{ven},i}^{\mathrm{sys}}, \quad\scriptstyle{i\in\{\mathrm{spl,espl,\atop msc,cer,cor}\}} && \text{systemic capillary arterial momentum}\nonumber\\ &C_{\mathrm{ven},i}^{\mathrm{sys}} \frac{\mathrm{d}p_{\mathrm{ven},i}^{\mathrm{sys}}}{\mathrm{d}t} = q_{\mathrm{ar},i}^{\mathrm{sys}} - q_{\mathrm{ven},i}^{\mathrm{sys}}, \quad\scriptstyle{i\in\{\mathrm{spl,espl,\atop msc,cer,cor}\}}&& \text{systemic capillary venous flow balance}\nonumber\\ &R_{\mathrm{ven},i}^{\mathrm{sys}}\,q_{\mathrm{ven},i}^{\mathrm{sys}} = p_{\mathrm{ven},i}^{\mathrm{sys}}-p_{\mathrm{ven}}^{\mathrm{sys}}, \quad\scriptstyle{i\in\{\mathrm{spl,espl,\atop msc,cer,cor}\}} && \text{systemic capillary venous momentum}\nonumber\\ &C_{\mathrm{ven}}^{\mathrm{sys}} \frac{\mathrm{d}p_{\mathrm{ven}}^{\mathrm{sys}}}{\mathrm{d}t} = \!\!\!\!\sum_{j=\mathrm{spl,espl,\atop msc,cer,cor}}\!\!\!\!\!q_{\mathrm{ven},j}^{\mathrm{sys}}-q_{\mathrm{ven}}^{\mathrm{sys}} && \text{systemic venous flow balance}\nonumber\\ &L_{\mathrm{ven}}^{\mathrm{sys}}\frac{\mathrm{d}q_{\mathrm{ven}}^{\mathrm{sys}}}{\mathrm{d}t} + R_{\mathrm{ven}}^{\mathrm{sys}}\, q_{\mathrm{ven}}^{\mathrm{sys}} = p_{\mathrm{ven}}^{\mathrm{sys}} - p_{\mathrm{at}}^{r} && \text{systemic venous momentum}\nonumber \end{aligned}\end{split}\]

(41)\[\begin{split}\begin{aligned} &-Q_{\mathrm{at}}^{r} = q_{\mathrm{ven}}^{\mathrm{sys}} - q_{\mathrm{v,in}}^{r} && \text{right atrium flow balance}\\ &\tilde{R}_{\mathrm{v,in}}^{r}\,q_{\mathrm{v,in}}^{r} = p_{\mathrm{at}}^{r}-p_{\mathrm{v}}^{r} && \text{tricuspid valve momentum}\\ &-Q_{\mathrm{v}}^{r} = q_{\mathrm{v,in}}^{r} - q_{\mathrm{v,out}}^{r} && \text{right ventricle flow balance}\nonumber\\ &\tilde{R}_{\mathrm{v,out}}^{r}\,q_{\mathrm{v,out}}^{r} = p_{\mathrm{v}}^{r}-p_{\mathrm{ar}}^{\mathrm{pul}} && \text{pulmonary valve momentum}\nonumber\\ &C_{\mathrm{ar}}^{\mathrm{pul}} \frac{\mathrm{d}p_{\mathrm{ar}}^{\mathrm{pul}}}{\mathrm{d}t} = q_{\mathrm{v,out}}^{r} - q_{\mathrm{ar}}^{\mathrm{pul}} && \text{pulmonary arterial flow balance}\\ &L_{\mathrm{ar}}^{\mathrm{pul}}\frac{\mathrm{d}q_{\mathrm{ar}}^{\mathrm{pul}}}{\mathrm{d}t} + R_{\mathrm{ar}}^{\mathrm{pul}}\,q_{\mathrm{ar}}^{\mathrm{pul}}=p_{\mathrm{ar}}^{\mathrm{pul}} -p_{\mathrm{cap}}^{\mathrm{pul}} && \text{pulmonary arterial momentum}\\ &C_{\mathrm{cap}}^{\mathrm{pul}} \frac{\mathrm{d}p_{\mathrm{cap}}^{\mathrm{pul}}}{\mathrm{d}t} = q_{\mathrm{ar}}^{\mathrm{pul}} - q_{\mathrm{cap}}^{\mathrm{pul}} && \text{pulmonary capillary flow balance}\\ &R_{\mathrm{cap}}^{\mathrm{pul}}\,q_{\mathrm{cap}}^{\mathrm{pul}}=p_{\mathrm{cap}}^{\mathrm{pul}}-p_{\mathrm{ven}}^{\mathrm{pul}} && \text{pulmonary capillary momentum}\\ &C_{\mathrm{ven}}^{\mathrm{pul}} \frac{\mathrm{d}p_{\mathrm{ven}}^{\mathrm{pul}}}{\mathrm{d}t} = q_{\mathrm{cap}}^{\mathrm{pul}} - q_{\mathrm{ven}}^{\mathrm{pul}} && \text{pulmonary venous flow balance}\\ &L_{\mathrm{ven}}^{\mathrm{pul}}\frac{\mathrm{d}q_{\mathrm{ven}}^{\mathrm{pul}}}{\mathrm{d}t} + R_{\mathrm{ven}}^{\mathrm{pul}}\, q_{\mathrm{ven}}^{\mathrm{pul}}=p_{\mathrm{ven}}^{\mathrm{pul}}-p_{\mathrm{at}}^{\ell} && \text{pulmonary venous momentum} \end{aligned}\end{split}\]

with:

\[\begin{aligned} Q_{\mathrm{at}}^{\ell} := -\frac{\mathrm{d}V_{\mathrm{at}}^{\ell}}{\mathrm{d}t}, \qquad Q_{\mathrm{v}}^{\ell} := -\frac{\mathrm{d}V_{\mathrm{v}}^{\ell}}{\mathrm{d}t}, \qquad Q_{\mathrm{at}}^{r} := -\frac{\mathrm{d}V_{\mathrm{at}}^{r}}{\mathrm{d}t}, \qquad Q_{\mathrm{v}}^{r} := -\frac{\mathrm{d}V_{\mathrm{v}}^{r}}{\mathrm{d}t}\nonumber \end{aligned}\]

Systemic and pulmonary circulation, including capillary and coronary flow

– Model type : syspulcapcor

– Variant of syspulcap, with coronaries branching off directly after aortic valve

(42)\[\begin{split}\begin{aligned} &-Q_{\mathrm{at}}^{\ell} = q_{\mathrm{ven}}^{\mathrm{pul}} - q_{\mathrm{v,in}}^{\ell} && \text{left atrium flow balance}\\ &\tilde{R}_{\mathrm{v,in}}^{\ell}\,q_{\mathrm{v,in}}^{\ell} = p_{\mathrm{at}}^{\ell}-p_{\mathrm{v}}^{\ell} && \text{mitral valve momentum}\nonumber\\ &-Q_{\mathrm{v}}^{\ell} = q_{\mathrm{v,in}}^{\ell} - q_{\mathrm{v,out}}^{\ell} && \text{left ventricle flow balance}\\ &\tilde{R}_{\mathrm{v,out}}^{\ell}\,q_{\mathrm{v,out}}^{\ell} = p_{\mathrm{v}}^{\ell}-p_{\mathrm{ar}}^{\mathrm{sys}} && \text{aortic valve momentum}\nonumber\\ &0 = q_{\mathrm{v,out}}^{\ell} - q_{\mathrm{ar,p}}^{\mathrm{sys}} - q_{\mathrm{ar,cor,in}}^{\mathrm{sys}} && \text{aortic root flow balance}\\ &I_{\mathrm{ar}}^{\mathrm{sys}} \frac{\mathrm{d}q_{\mathrm{ar,p}}^{\mathrm{sys}}}{\mathrm{d}t} + Z_{\mathrm{ar}}^{\mathrm{sys}}\,q_{\mathrm{ar,p}}^{\mathrm{sys}}=p_{\mathrm{ar}}^{\mathrm{sys}}-p_{\mathrm{ar,d}}^{\mathrm{sys}} && \text{aortic root inertia}\nonumber\\ &C_{\mathrm{ar,cor}}^{\mathrm{sys}} \frac{\mathrm{d}p_{\mathrm{ar}}^{\mathrm{sys}}}{\mathrm{d}t} = q_{\mathrm{ar,cor,in}}^{\mathrm{sys}} - q_{\mathrm{ar,cor}}^{\mathrm{sys}} && \text{systemic arterial coronary flow balance}\nonumber\\ &R_{\mathrm{ar,cor}}^{\mathrm{sys}}\,q_{\mathrm{ar,cor}}^{\mathrm{sys}} = p_{\mathrm{ar}}^{\mathrm{sys}} - p_{\mathrm{ven,cor}}^{\mathrm{sys}} && \text{systemic arterial coronary momentum}\\ &C_{\mathrm{ar}}^{\mathrm{sys}} \frac{\mathrm{d}p_{\mathrm{ar,d}}^{\mathrm{sys}}}{\mathrm{d}t} = q_{\mathrm{ar,p}}^{\mathrm{sys}} - q_{\mathrm{ar}}^{\mathrm{sys}}&& \text{systemic arterial flow balance}\nonumber\\ &L_{\mathrm{ar}}^{\mathrm{sys}}\frac{\mathrm{d}q_{\mathrm{ar}}^{\mathrm{sys}}}{\mathrm{d}t} + R_{\mathrm{ar}}^{\mathrm{sys}}\,q_{\mathrm{ar}}^{\mathrm{sys}}=p_{\mathrm{ar,d}}^{\mathrm{sys}} -p_{\mathrm{ar,peri}}^{\mathrm{sys}} && \text{systemic arterial flow balance}\\ &\left(\sum_{j\in\{\mathrm{spl,espl,\atop msc,cer}\}}\!\!\!\!\!\!\!\!\!C_{\mathrm{ar},j}^{\mathrm{sys}}\right) \frac{\mathrm{d}p_{\mathrm{ar,peri}}^{\mathrm{sys}}}{\mathrm{d}t} = q_{\mathrm{ar}}^{\mathrm{sys}}-\!\!\!\!\!\sum_{j\in\{\mathrm{spl,espl,\atop msc,cer}\}}\!\!\!\!\!\!\!\!\!q_{\mathrm{ar},j}^{\mathrm{sys}} && \text{systemic arterial capillary flow balance}\\ &R_{\mathrm{ar},i}^{\mathrm{sys}}\,q_{\mathrm{ar},i}^{\mathrm{sys}} = p_{\mathrm{ar,peri}}^{\mathrm{sys}} - p_{\mathrm{ven},i}^{\mathrm{sys}}, \quad\scriptstyle{i\in\{\mathrm{spl,espl,\atop msc,cer}\}} && \text{systemic arterial capillary momentum}\\ &C_{\mathrm{ven},i}^{\mathrm{sys}} \frac{\mathrm{d}p_{\mathrm{ven},i}^{\mathrm{sys}}}{\mathrm{d}t} = q_{\mathrm{ar},i}^{\mathrm{sys}} - q_{\mathrm{ven},i}^{\mathrm{sys}}, \quad\scriptstyle{i\in\{\mathrm{spl,espl,\atop msc,cer}\}} && \text{systemic venous capillary flow balance}\nonumber\\ &R_{\mathrm{ven},i}^{\mathrm{sys}}\,q_{\mathrm{ven},i}^{\mathrm{sys}} = p_{\mathrm{ven},i}^{\mathrm{sys}}-p_{\mathrm{ven}}^{\mathrm{sys}}, \quad\scriptstyle{i\in\{\mathrm{spl,espl,\atop msc,cer}\}} && \text{systemic venous capillary momentum}\nonumber\\ &C_{\mathrm{ven}}^{\mathrm{sys}} \frac{\mathrm{d}p_{\mathrm{ven}}^{\mathrm{sys}}}{\mathrm{d}t} = \!\!\!\!\sum_{j=\mathrm{spl,espl,\atop msc,cer}}\!\!\!\!\!q_{\mathrm{ven},j}^{\mathrm{sys}}-q_{\mathrm{ven}}^{\mathrm{sys}} && \text{systemic venous flow balance}\\ &L_{\mathrm{ven}}^{\mathrm{sys}}\frac{\mathrm{d}q_{\mathrm{ven}}^{\mathrm{sys}}}{\mathrm{d}t} + R_{\mathrm{ven}}^{\mathrm{sys}}\, q_{\mathrm{ven}}^{\mathrm{sys}} = p_{\mathrm{ven}}^{\mathrm{sys}} - p_{\mathrm{at}}^{r} && \text{systemic venous momentum}\nonumber\\ &C_{\mathrm{ven,cor}}^{\mathrm{sys}} \frac{\mathrm{d}p_{\mathrm{ven,cor}}^{\mathrm{sys}}}{\mathrm{d}t} = q_{\mathrm{ar,cor}}^{\mathrm{sys}}-q_{\mathrm{ven,cor}}^{\mathrm{sys}} && \text{systemic venous coronary flow balance}\nonumber\\ &R_{\mathrm{ven,cor}}^{\mathrm{sys}}\,q_{\mathrm{ven,cor}}^{\mathrm{sys}} = p_{\mathrm{ven,cor}}^{\mathrm{sys}} - p_{\mathrm{at}}^{r} && \text{systemic venous coronary momentum} \end{aligned}\end{split}\]

(43)\[\begin{split}\begin{aligned} &-Q_{\mathrm{at}}^{r} = q_{\mathrm{ven}}^{\mathrm{sys}} + q_{\mathrm{ven,cor}}^{\mathrm{sys}} - q_{\mathrm{v,in}}^{r} && \text{right atrium flow balance}\nonumber\\ &\tilde{R}_{\mathrm{v,in}}^{r}\,q_{\mathrm{v,in}}^{r} = p_{\mathrm{at}}^{r}-p_{\mathrm{v}}^{r} && \text{tricuspid valve momentum}\\ &-Q_{\mathrm{v}}^{r} = q_{\mathrm{v,in}}^{r} - q_{\mathrm{v,out}}^{r} && \text{right ventricle flow balance}\\ &\tilde{R}_{\mathrm{v,out}}^{r}\,q_{\mathrm{v,out}}^{r} = p_{\mathrm{v}}^{r}-p_{\mathrm{ar}}^{\mathrm{pul}} && \text{pulmonary valve momentum}\\ &C_{\mathrm{ar}}^{\mathrm{pul}} \frac{\mathrm{d}p_{\mathrm{ar}}^{\mathrm{pul}}}{\mathrm{d}t} = q_{\mathrm{v,out}}^{r} - q_{\mathrm{ar}}^{\mathrm{pul}} && \text{pulmonary arterial flow balance}\\ &L_{\mathrm{ar}}^{\mathrm{pul}}\frac{\mathrm{d}q_{\mathrm{ar}}^{\mathrm{pul}}}{\mathrm{d}t} + R_{\mathrm{ar}}^{\mathrm{pul}}\,q_{\mathrm{ar}}^{\mathrm{pul}}=p_{\mathrm{ar}}^{\mathrm{pul}} -p_{\mathrm{cap}}^{\mathrm{pul}} && \text{pulmonary arterial momentum}\\ &C_{\mathrm{cap}}^{\mathrm{pul}} \frac{\mathrm{d}p_{\mathrm{cap}}^{\mathrm{pul}}}{\mathrm{d}t} = q_{\mathrm{ar}}^{\mathrm{pul}} - q_{\mathrm{cap}}^{\mathrm{pul}} && \text{pulmonary capillary flow balance}\\ &R_{\mathrm{cap}}^{\mathrm{pul}}\,q_{\mathrm{cap}}^{\mathrm{pul}}=p_{\mathrm{cap}}^{\mathrm{pul}}-p_{\mathrm{ven}}^{\mathrm{pul}} && \text{pulmonary capillary momentum}\\ &C_{\mathrm{ven}}^{\mathrm{pul}} \frac{\mathrm{d}p_{\mathrm{ven}}^{\mathrm{pul}}}{\mathrm{d}t} = q_{\mathrm{cap}}^{\mathrm{pul}} - q_{\mathrm{ven}}^{\mathrm{pul}} && \text{pulmonary venous flow balance}\\ &L_{\mathrm{ven}}^{\mathrm{pul}}\frac{\mathrm{d}q_{\mathrm{ven}}^{\mathrm{pul}}}{\mathrm{d}t} + R_{\mathrm{ven}}^{\mathrm{pul}}\, q_{\mathrm{ven}}^{\mathrm{pul}}=p_{\mathrm{ven}}^{\mathrm{pul}}-p_{\mathrm{at}}^{\ell} && \text{pulmonary venous momentum} \end{aligned}\end{split}\]

with:

\[\begin{aligned} Q_{\mathrm{at}}^{\ell} := -\frac{\mathrm{d}V_{\mathrm{at}}^{\ell}}{\mathrm{d}t}, \qquad Q_{\mathrm{v}}^{\ell} := -\frac{\mathrm{d}V_{\mathrm{v}}^{\ell}}{\mathrm{d}t}, \qquad Q_{\mathrm{at}}^{r} := -\frac{\mathrm{d}V_{\mathrm{at}}^{r}}{\mathrm{d}t}, \qquad Q_{\mathrm{v}}^{r} := -\frac{\mathrm{d}V_{\mathrm{v}}^{r}}{\mathrm{d}t} \end{aligned}\]

Systemic and pulmonary circulation, capillary flow, and respiratory (gas transport + dissociation) model

– Model type : syspulcaprespir

– Model equations described in [Hir19], p. 51ff., 58ff.

Multi-physics coupling

Solid + 0D flow

– Example: Sec. 4.4.1 and demos/solid_flow0d

– Problem type: solid_flow0d

– solid momentum in ([equation-solid-weak-form]) or in ([equation-solid-weak-form-inc]) augmented by following term:

\[\begin{aligned} r_u \leftarrow r_u + \int\limits_{\mathit{\Gamma}_0^{\text{s}\text{-}\mathrm{0d}}}\!\mathit{\Lambda}\,J\boldsymbol{F}^{-\mathrm{T}}\boldsymbol{n}_0\cdot\delta\boldsymbol{u}\,\mathrm{d}A \end{aligned}\]

– Multiplier constraint

\[\begin{aligned} r_{\mathit{\Lambda}}(\mathit{\Lambda},\boldsymbol{u};\delta\mathit{\Lambda}):= \left(\int\limits_{\mathit{\Gamma}_0^{\mathrm{\text{s}\text{-}0d}}}\! J\boldsymbol{F}^{-\mathrm{T}}\boldsymbol{n}_{0}\cdot\boldsymbol{v}(\boldsymbol{u})\,\mathrm{d}A - Q^{\mathrm{0d}}(\mathit{\Lambda})\right) \delta\mathit{\Lambda}, \quad \forall \; \delta\mathit{\Lambda} \end{aligned}\]

– Discrete nonlinear system to solve in each time step \(n\) for displacement-based solid:

(44)\[\begin{split}\begin{aligned} \boldsymbol{\mathsf{r}}_{n+1} = \begin{bmatrix} \boldsymbol{\mathsf{r}}_{u}(\boldsymbol{\mathsf{u}},\boldsymbol{\mathsf{\Lambda}}) \\ \boldsymbol{\mathsf{r}}_{\mathit{\Lambda}}(\boldsymbol{\mathsf{\Lambda}},\boldsymbol{\mathsf{u}}) \end{bmatrix}_{n+1} = \boldsymbol{\mathsf{0}} \end{aligned}\end{split}\]

– Discrete linear system to solve in each Newton iteration \(k\) for displacement-based solid:

(45)\[\begin{split}\begin{aligned} \begin{bmatrix} \boldsymbol{\mathsf{K}}_{uu} & \boldsymbol{\mathsf{K}}_{u\mathit{\Lambda}} \\ \\ \boldsymbol{\mathsf{K}}_{\mathit{\Lambda}u} & \boldsymbol{\mathsf{K}}_{\mathit{\Lambda}\mathit{\Lambda}}\end{bmatrix}_{n+1}^{k}\begin{bmatrix} \Delta\boldsymbol{\mathsf{u}} \\ \\ \Delta\boldsymbol{\mathsf{\Lambda}}\end{bmatrix}_{n+1}^{k+1}=-\begin{bmatrix} \boldsymbol{\mathsf{r}}_{u} \\ \\ \boldsymbol{\mathsf{r}}_{\mathit{\Lambda}}\end{bmatrix}_{n+1}^{k} \end{aligned}\end{split}\]

– Discrete nonlinear system to solve in each time step \(n\) for incompressible solid:

(46)\[\begin{split}\begin{aligned} \boldsymbol{\mathsf{r}}_{n+1} = \begin{bmatrix} \boldsymbol{\mathsf{r}}_{u}(\boldsymbol{\mathsf{u}},\boldsymbol{\mathsf{p}},\boldsymbol{\mathsf{\Lambda}}) \\ \boldsymbol{\mathsf{r}}_{p}(\boldsymbol{\mathsf{u}}) \\ \boldsymbol{\mathsf{r}}_{\mathit{\Lambda}}(\boldsymbol{\mathsf{\Lambda}},\boldsymbol{\mathsf{u}}) \end{bmatrix}_{n+1} = \boldsymbol{\mathsf{0}} \end{aligned}\end{split}\]

– Discrete linear system to solve in each Newton iteration \(k\) for incompressible solid:

(47)\[\begin{split}\begin{aligned} \begin{bmatrix} \boldsymbol{\mathsf{K}}_{uu} & \boldsymbol{\mathsf{K}}_{up} & \boldsymbol{\mathsf{K}}_{u\mathit{\Lambda}} \\ \\ \boldsymbol{\mathsf{K}}_{pu} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}\\ \\ \boldsymbol{\mathsf{K}}_{\mathit{\Lambda}u} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{\mathit{\Lambda}\mathit{\Lambda}}\end{bmatrix}_{n+1}^{k}\begin{bmatrix} \Delta\boldsymbol{\mathsf{u}} \\ \\ \Delta\boldsymbol{\mathsf{p}} \\ \\ \Delta\boldsymbol{\mathsf{\Lambda}}\end{bmatrix}_{n+1}^{k+1}=-\begin{bmatrix} \boldsymbol{\mathsf{r}}_{u} \\ \\ \boldsymbol{\mathsf{r}}_{p} \\ \\ \boldsymbol{\mathsf{r}}_{\mathit{\Lambda}}\end{bmatrix}_{n+1}^{k} \end{aligned}\end{split}\]

– sub-solves: 0D model has to hold true in each nonlinear iteration \(k\):

(48)\[\begin{aligned} \left.\boldsymbol{\mathsf{r}}^{0\mathrm{D}}(\boldsymbol{\mathsf{\Lambda}})\right|_{k} = \boldsymbol{\mathsf{0}}, \qquad \rightsquigarrow \boldsymbol{\mathsf{K}}_{k}^{0\mathrm{D},j}\Delta\boldsymbol{\mathsf{q}}_{k}^{j+1} = -\boldsymbol{\mathsf{r}}_{k}^{0\mathrm{D},j} \end{aligned}\]

Fluid + 0D flow

– Example: Sec. 4.4.2 demos/fluid_flow0d

– Problem type: fluid_flow0d

– fluid momentum in ([equation-fluid-weak-form]) augmented by following term:

\[\begin{aligned} r_v \leftarrow r_v + \int\limits_{\mathit{\Gamma}_t^{\text{f}\text{-}\mathrm{0d}}}\!\mathit{\Lambda}\,\boldsymbol{n}\cdot\delta\boldsymbol{v}\,\mathrm{d}a \end{aligned}\]

– Multiplier constraint

\[\begin{aligned} r_{\mathit{\Lambda}}(\mathit{\Lambda},\boldsymbol{v};\delta\mathit{\Lambda}):= \left(\int\limits_{\mathit{\Gamma}_t^{\mathrm{\text{f}\text{-}0d}}}\! \boldsymbol{n}\cdot\boldsymbol{v}\,\mathrm{d}a - Q^{\mathrm{0d}}(\mathit{\Lambda})\right) \delta\mathit{\Lambda}, \quad \forall \; \delta\mathit{\Lambda} \end{aligned}\]

– Discrete nonlinear system to solve in each time step \(n\):

(49)\[\begin{split}\begin{aligned} \boldsymbol{\mathsf{r}}_{n+1} = \begin{bmatrix} \boldsymbol{\mathsf{r}}_{v}(\boldsymbol{\mathsf{v}},\boldsymbol{\mathsf{p}},\boldsymbol{\mathsf{\Lambda}}) \\ \boldsymbol{\mathsf{r}}_{p}(\boldsymbol{\mathsf{p}},\boldsymbol{\mathsf{v}}) \\ \boldsymbol{\mathsf{r}}_{\mathit{\Lambda}}(\boldsymbol{\mathsf{\Lambda}},\boldsymbol{\mathsf{v}}) \end{bmatrix}_{n+1} = \boldsymbol{\mathsf{0}} \end{aligned}\end{split}\]

– Discrete linear system to solve in each Newton iteration \(k\):

(50)\[\begin{split}\begin{aligned} \begin{bmatrix} \boldsymbol{\mathsf{K}}_{vv} & \boldsymbol{\mathsf{K}}_{vp} & \boldsymbol{\mathsf{K}}_{v\mathit{\Lambda}} \\ \\ \boldsymbol{\mathsf{K}}_{pv} & \boldsymbol{\mathsf{K}}_{pp} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}\\ \\ \boldsymbol{\mathsf{K}}_{\mathit{\Lambda}v} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{\mathit{\Lambda}\mathit{\Lambda}}\end{bmatrix}_{n+1}^{k}\begin{bmatrix} \Delta\boldsymbol{\mathsf{v}} \\ \\ \Delta\boldsymbol{\mathsf{p}} \\ \\ \Delta\boldsymbol{\mathsf{\Lambda}}\end{bmatrix}_{n+1}^{k+1}=-\begin{bmatrix} \boldsymbol{\mathsf{r}}_{v} \\ \\ \boldsymbol{\mathsf{r}}_{p} \\ \\ \boldsymbol{\mathsf{r}}_{\mathit{\Lambda}}\end{bmatrix}_{n+1}^{k} \end{aligned}\end{split}\]

ALE fluid + 0D flow

– Problem type: fluid_ale_flow0d

– fluid momentum in ([equation-fluid-ale-weak-form]) augmented by following term:

\[\begin{aligned} r_v \leftarrow r_v + \int\limits_{\mathit{\Gamma}_0^{\text{f}\text{-}\mathrm{0d}}}\!\mathit{\Lambda}\,\widehat{J}\widehat{\boldsymbol{F}}^{-\mathrm{T}}\boldsymbol{n}_{0}\cdot\delta\boldsymbol{v}\,\mathrm{d}A \end{aligned}\]

– Multiplier constraint

\[\begin{aligned} r_{\mathit{\Lambda}}(\mathit{\Lambda},\boldsymbol{v},\boldsymbol{d};\delta\mathit{\Lambda}):= \left(\int\limits_{\mathit{\Gamma}_0^{\mathrm{\text{f}\text{-}0d}}}\! \widehat{J}\widehat{\boldsymbol{F}}^{-\mathrm{T}}\boldsymbol{n}_{0}\cdot(\boldsymbol{v}-\widehat{\boldsymbol{w}}(\boldsymbol{d}))\,\mathrm{d}A - Q^{\mathrm{0d}}(\mathit{\Lambda})\right) \delta\mathit{\Lambda}, \quad \forall \; \delta\mathit{\Lambda} \end{aligned}\]

with \(\boldsymbol{w}(\boldsymbol{d})=\frac{\mathrm{d}\boldsymbol{d}}{\mathrm{d}t}\)

– Discrete nonlinear system to solve in each time step \(n\):

(51)\[\begin{split}\begin{aligned} \boldsymbol{\mathsf{r}}_{n+1} = \begin{bmatrix} \boldsymbol{\mathsf{r}}_{v}(\boldsymbol{\mathsf{v}},\boldsymbol{\mathsf{p}},\boldsymbol{\mathsf{\Lambda}},\boldsymbol{\mathsf{d}}) \\ \boldsymbol{\mathsf{r}}_{p}(\boldsymbol{\mathsf{p}},\boldsymbol{\mathsf{v}},\boldsymbol{\mathsf{d}}) \\ \boldsymbol{\mathsf{r}}_{\mathit{\Lambda}}(\boldsymbol{\mathsf{\Lambda}},\boldsymbol{\mathsf{v}},\boldsymbol{\mathsf{d}}) \\ \boldsymbol{\mathsf{r}}_{d}(\boldsymbol{\mathsf{d}}) \end{bmatrix}_{n+1} = \boldsymbol{\mathsf{0}} \end{aligned}\end{split}\]

– Discrete linear system to solve in each Newton iteration \(k\):

(52)\[\begin{split}\begin{aligned} \begin{bmatrix} \boldsymbol{\mathsf{K}}_{vv} & \boldsymbol{\mathsf{K}}_{vp} & \boldsymbol{\mathsf{K}}_{v\mathit{\Lambda}} & \boldsymbol{\mathsf{K}}_{vd} \\ \\ \boldsymbol{\mathsf{K}}_{pv} & \boldsymbol{\mathsf{K}}_{pp} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{pd} \\ \\ \boldsymbol{\mathsf{K}}_{\mathit{\Lambda}v} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{\mathit{\Lambda}\mathit{\Lambda}} & \boldsymbol{\mathsf{K}}_{\mathit{\Lambda}d} \\ \\ \boldsymbol{\mathsf{K}}_{dv} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{dd} \end{bmatrix}_{n+1}^{k}\begin{bmatrix} \Delta\boldsymbol{\mathsf{v}} \\ \\ \Delta\boldsymbol{\mathsf{p}} \\ \\ \Delta\boldsymbol{\mathsf{\Lambda}}\\ \\ \Delta\boldsymbol{\mathsf{d}} \end{bmatrix}_{n+1}^{k+1}=-\begin{bmatrix} \boldsymbol{\mathsf{r}}_{v} \\ \\ \boldsymbol{\mathsf{r}}_{p} \\ \\ \boldsymbol{\mathsf{r}}_{\mathit{\Lambda}} \\ \\ \boldsymbol{\mathsf{r}}_{d}\end{bmatrix}_{n+1}^{k} \end{aligned}\end{split}\]

Fluid-Solid Interaction (FSI)

– Example: Sec. 5.6 and demos/fsi

– Problem type: fsi

– Solid momentum in ([equation-solid-weak-form]) or ([equation-solid-weak-form-inc]) augmented by following term:

\[\begin{aligned} r_u \leftarrow r_u + \int\limits_{\mathit{\Gamma}_0^{\text{f}\text{-}\text{s}}}\boldsymbol{\lambda}\cdot\delta\boldsymbol{u}\,\mathrm{d}A \end{aligned}\]

– Fluid momentum in ([equation-fluid-ale-weak-form]) is augmented by following term:

\[\begin{aligned} r_v \leftarrow r_v - \int\limits_{\mathit{\Gamma}_0^{\text{f}\text{-}\text{s}}}\boldsymbol{\lambda}\cdot\delta\boldsymbol{v}\,\mathrm{d}A \end{aligned}\]

Note the different signs (actio=reactio!)

– Lagrange multiplier constraint:

• “solid-governed”:

  \[\begin{aligned} r_{\lambda}(\boldsymbol{u},\boldsymbol{v};\delta\boldsymbol{\lambda}):= \int\limits_{\mathit{\Gamma}_0^{\mathrm{\text{f}\text{-}\text{s}}}} \left(\boldsymbol{u} - \boldsymbol{u}_{\mathrm{f}}(\boldsymbol{v})\right)\cdot\delta\boldsymbol{\lambda}\,\mathrm{d}A, \quad \forall \; \delta\boldsymbol{\lambda} \end{aligned}\]

• “fluid-governed”:

  \[\begin{aligned} r_{\lambda}(\boldsymbol{v},\boldsymbol{u};\delta\boldsymbol{\lambda}):= \int\limits_{\mathit{\Gamma}_0^{\mathrm{\text{f}\text{-}\text{s}}}} \left(\boldsymbol{v} - \frac{\mathrm{d} \boldsymbol{u}}{\mathrm{d} t}\right)\cdot\delta\boldsymbol{\lambda}\,\mathrm{d}A, \quad \forall \; \delta\boldsymbol{\lambda} \end{aligned}\]

– Discrete nonlinear system to solve in each time step \(n\) for displacement-based solid:

(53)\[\begin{split}\begin{aligned} \boldsymbol{\mathsf{r}}_{n+1} = \begin{bmatrix} \boldsymbol{\mathsf{r}}_{u}^{\mathrm{s}}(\boldsymbol{\mathsf{u}},\boldsymbol{\mathsf{\lambda}}) \\ \boldsymbol{\mathsf{r}}_{v}^{\mathrm{f}}(\boldsymbol{\mathsf{v}},\boldsymbol{\mathsf{p}},\boldsymbol{\mathsf{\lambda}},\boldsymbol{\mathsf{d}}) \\ \boldsymbol{\mathsf{r}}_{p}^{\mathrm{f}}(\boldsymbol{\mathsf{p}},\boldsymbol{\mathsf{v}},\boldsymbol{\mathsf{d}}) \\ \boldsymbol{\mathsf{r}}_{\lambda}(\boldsymbol{\mathsf{u}},\boldsymbol{\mathsf{v}}) \\ \boldsymbol{\mathsf{r}}_{d}(\boldsymbol{\mathsf{d}}) \end{bmatrix}_{n+1} = \boldsymbol{\mathsf{0}} \end{aligned}\end{split}\]

– Discrete nonlinear system to solve in each time step \(n\) for incompressible solid:

(54)\[\begin{split}\begin{aligned} \boldsymbol{\mathsf{r}}_{n+1} = \begin{bmatrix} \boldsymbol{\mathsf{r}}_{u}^{\mathrm{s}}(\boldsymbol{\mathsf{u}},\boldsymbol{\mathsf{p}}^{\mathrm{s}},\boldsymbol{\mathsf{\lambda}}) \\ \boldsymbol{\mathsf{r}}_{p}^{\mathrm{s}}(\boldsymbol{\mathsf{u}}) \\ \boldsymbol{\mathsf{r}}_{v}^{\mathrm{f}}(\boldsymbol{\mathsf{v}},\boldsymbol{\mathsf{p}},\boldsymbol{\mathsf{\lambda}},\boldsymbol{\mathsf{d}}) \\ \boldsymbol{\mathsf{r}}_{p}^{\mathrm{f}}(\boldsymbol{\mathsf{p}},\boldsymbol{\mathsf{v}},\boldsymbol{\mathsf{d}}) \\ \boldsymbol{\mathsf{r}}_{\lambda}(\boldsymbol{\mathsf{u}},\boldsymbol{\mathsf{v}}) \\ \boldsymbol{\mathsf{r}}_{d}(\boldsymbol{\mathsf{d}}) \end{bmatrix}_{n+1} = \boldsymbol{\mathsf{0}} \end{aligned}\end{split}\]

– Discrete linear system to solve in each Newton iteration \(k\) for displacement-based solid:

(55)\[\begin{split}\begin{aligned} \begin{bmatrix} \boldsymbol{\mathsf{K}}_{uu} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{u\lambda} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}\\ \\ \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{vv} & \boldsymbol{\mathsf{K}}_{vp} & \boldsymbol{\mathsf{K}}_{v\lambda} & \boldsymbol{\mathsf{K}}_{vd} \\ \\ \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{pv} & \boldsymbol{\mathsf{K}}_{pp} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{pd} \\ \\ \boldsymbol{\mathsf{K}}_{\lambda u} & \boldsymbol{\mathsf{K}}_{\lambda v} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}\\ \\ \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{dv} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{dd} \end{bmatrix}_{n+1}^{k}\begin{bmatrix} \Delta\boldsymbol{\mathsf{u}} \\ \\ \Delta\boldsymbol{\mathsf{v}} \\ \\ \Delta\boldsymbol{\mathsf{p}} \\ \\ \Delta\boldsymbol{\mathsf{\lambda}}\\ \\ \Delta\boldsymbol{\mathsf{d}} \end{bmatrix}_{n+1}^{k+1}=-\begin{bmatrix} \boldsymbol{\mathsf{r}}_{u}^{\mathrm{s}} \\ \\ \boldsymbol{\mathsf{r}}_{v}^{\mathrm{f}} \\ \\ \boldsymbol{\mathsf{r}}_{p}^{\mathrm{f}} \\ \\ \boldsymbol{\mathsf{r}}_{\lambda} \\ \\ \boldsymbol{\mathsf{r}}_{d}\end{bmatrix}_{n+1}^{k} \end{aligned}\end{split}\]

– Discrete linear system to solve in each Newton iteration \(k\) for incompressible solid:

(56)\[\begin{split}\begin{aligned} \begin{bmatrix} \boldsymbol{\mathsf{K}}_{uu} & \boldsymbol{\mathsf{K}}_{up} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{u\lambda} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}\\ \\ \boldsymbol{\mathsf{K}}_{pu} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}\\ \\ \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{vv} & \boldsymbol{\mathsf{K}}_{vp} & \boldsymbol{\mathsf{K}}_{v\lambda} & \boldsymbol{\mathsf{K}}_{vd} \\ \\ \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{pv} & \boldsymbol{\mathsf{K}}_{pp} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{pd} \\ \\ \boldsymbol{\mathsf{K}}_{\lambda u} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{\lambda v} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}\\ \\ \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{dv} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{dd} \end{bmatrix}_{n+1}^{k}\begin{bmatrix} \Delta\boldsymbol{\mathsf{u}} \\ \\ \Delta\boldsymbol{\mathsf{p}} \\ \\ \Delta\boldsymbol{\mathsf{v}} \\ \\ \Delta\boldsymbol{\mathsf{p}} \\ \\ \Delta\boldsymbol{\mathsf{\lambda}}\\ \\ \Delta\boldsymbol{\mathsf{d}} \end{bmatrix}_{n+1}^{k+1}=-\begin{bmatrix} \boldsymbol{\mathsf{r}}_{u}^{\mathrm{s}} \\ \\ \boldsymbol{\mathsf{r}}_{p}^{\mathrm{s}} \\ \\ \boldsymbol{\mathsf{r}}_{v}^{\mathrm{f}} \\ \\ \boldsymbol{\mathsf{r}}_{p}^{\mathrm{f}} \\ \\ \boldsymbol{\mathsf{r}}_{\lambda} \\ \\ \boldsymbol{\mathsf{r}}_{d}\end{bmatrix}_{n+1}^{k} \end{aligned}\end{split}\]

Fluid-Solid Interaction (FSI) + 0D flow

– Problem type: fsi_flow0d

– Discrete nonlinear system to solve in each time step \(n\) for displacement-based solid:

(57)\[\begin{split}\begin{aligned} \boldsymbol{\mathsf{r}}_{n+1} = \begin{bmatrix} \boldsymbol{\mathsf{r}}_{u}^{\mathrm{s}}(\boldsymbol{\mathsf{u}},\boldsymbol{\mathsf{\lambda}}) \\ \boldsymbol{\mathsf{r}}_{v}^{\mathrm{f}}(\boldsymbol{\mathsf{v}},\boldsymbol{\mathsf{p}},\boldsymbol{\mathsf{\lambda}},\boldsymbol{\mathsf{\Lambda}},\boldsymbol{\mathsf{d}}) \\ \boldsymbol{\mathsf{r}}_{p}^{\mathrm{f}}(\boldsymbol{\mathsf{p}},\boldsymbol{\mathsf{v}},\boldsymbol{\mathsf{d}}) \\ \boldsymbol{\mathsf{r}}_{\lambda}(\boldsymbol{\mathsf{u}},\boldsymbol{\mathsf{v}}) \\ \boldsymbol{\mathsf{r}}_{\mathit{\Lambda}}(\boldsymbol{\mathsf{\Lambda}},\boldsymbol{\mathsf{v}},\boldsymbol{\mathsf{d}}) \\ \boldsymbol{\mathsf{r}}_{d}(\boldsymbol{\mathsf{d}}) \end{bmatrix}_{n+1} = \boldsymbol{\mathsf{0}} \end{aligned}\end{split}\]

– Discrete nonlinear system to solve in each time step \(n\) for incompressible solid:

(58)\[\begin{split}\begin{aligned} \boldsymbol{\mathsf{r}}_{n+1} = \begin{bmatrix} \boldsymbol{\mathsf{r}}_{u}^{\mathrm{s}}(\boldsymbol{\mathsf{u}},\boldsymbol{\mathsf{p}}^{\mathrm{s}},\boldsymbol{\mathsf{\lambda}}) \\ \boldsymbol{\mathsf{r}}_{p}^{\mathrm{s}}(\boldsymbol{\mathsf{u}}) \\ \boldsymbol{\mathsf{r}}_{v}^{\mathrm{f}}(\boldsymbol{\mathsf{v}},\boldsymbol{\mathsf{p}},\boldsymbol{\mathsf{\lambda}},\boldsymbol{\mathsf{\Lambda}},\boldsymbol{\mathsf{d}}) \\ \boldsymbol{\mathsf{r}}_{p}^{\mathrm{f}}(\boldsymbol{\mathsf{p}},\boldsymbol{\mathsf{v}},\boldsymbol{\mathsf{d}}) \\ \boldsymbol{\mathsf{r}}_{\lambda}(\boldsymbol{\mathsf{u}},\boldsymbol{\mathsf{v}}) \\ \boldsymbol{\mathsf{r}}_{\mathit{\Lambda}}(\boldsymbol{\mathsf{\Lambda}},\boldsymbol{\mathsf{v}},\boldsymbol{\mathsf{d}}) \\ \boldsymbol{\mathsf{r}}_{d}(\boldsymbol{\mathsf{d}}) \end{bmatrix}_{n+1} = \boldsymbol{\mathsf{0}} \end{aligned}\end{split}\]

– Discrete linear system to solve in each Newton iteration \(k\) for displacement-based solid:

(59)\[\begin{split}\begin{aligned} \begin{bmatrix} \boldsymbol{\mathsf{K}}_{uu} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{u\lambda} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}\\ \\ \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{vv} & \boldsymbol{\mathsf{K}}_{vp} & \boldsymbol{\mathsf{K}}_{v\lambda} & \boldsymbol{\mathsf{K}}_{v\mathit{\Lambda}} & \boldsymbol{\mathsf{K}}_{vd} \\ \\ \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{pv} & \boldsymbol{\mathsf{K}}_{pp} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{pd} \\ \\ \boldsymbol{\mathsf{K}}_{\lambda u} & \boldsymbol{\mathsf{K}}_{\lambda v} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}\\ \\ \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{\mathit{\Lambda}v} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{\mathit{\Lambda}\mathit{\Lambda}} & \boldsymbol{\mathsf{K}}_{\mathit{\Lambda}d} \\ \\ \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{dv} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{dd} \end{bmatrix}_{n+1}^{k}\begin{bmatrix} \Delta\boldsymbol{\mathsf{u}} \\ \\ \Delta\boldsymbol{\mathsf{v}} \\ \\ \Delta\boldsymbol{\mathsf{p}} \\ \\ \Delta\boldsymbol{\mathsf{\lambda}}\\ \\ \Delta\boldsymbol{\mathsf{\Lambda}}\\ \\ \Delta\boldsymbol{\mathsf{d}} \end{bmatrix}_{n+1}^{k+1}=-\begin{bmatrix} \boldsymbol{\mathsf{r}}_{u}^{\mathrm{s}} \\ \\ \boldsymbol{\mathsf{r}}_{v}^{\mathrm{f}} \\ \\ \boldsymbol{\mathsf{r}}_{p}^{\mathrm{f}} \\ \\ \boldsymbol{\mathsf{r}}_{\lambda} \\ \\ \boldsymbol{\mathsf{r}}_{\mathit{\Lambda}} \\ \\ \boldsymbol{\mathsf{r}}_{d}\end{bmatrix}_{n+1}^{k} \end{aligned}\end{split}\]

– Discrete linear system to solve in each Newton iteration \(k\) for incompressible solid:

(60)\[\begin{split}\begin{aligned} \begin{bmatrix} \boldsymbol{\mathsf{K}}_{uu} & \boldsymbol{\mathsf{K}}_{up} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{u\lambda} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}\\ \\ \boldsymbol{\mathsf{K}}_{pu} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}\\ \\ \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{vv} & \boldsymbol{\mathsf{K}}_{vp} & \boldsymbol{\mathsf{K}}_{v\lambda} & \boldsymbol{\mathsf{K}}_{v\mathit{\Lambda}} & \boldsymbol{\mathsf{K}}_{vd} \\ \\ \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{pv} & \boldsymbol{\mathsf{K}}_{pp} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{pd} \\ \\ \boldsymbol{\mathsf{K}}_{\lambda u} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{\lambda v} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}\\ \\ \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{\mathit{\Lambda}v} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{\mathit{\Lambda}\mathit{\Lambda}} & \boldsymbol{\mathsf{K}}_{\mathit{\Lambda}d} \\ \\ \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{dv} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{dd} \end{bmatrix}_{n+1}^{k}\begin{bmatrix} \Delta\boldsymbol{\mathsf{u}} \\ \\ \Delta\boldsymbol{\mathsf{p}} \\ \\ \Delta\boldsymbol{\mathsf{v}} \\ \\ \Delta\boldsymbol{\mathsf{p}} \\ \\ \Delta\boldsymbol{\mathsf{\lambda}}\\ \\ \Delta\boldsymbol{\mathsf{\Lambda}}\\ \\ \Delta\boldsymbol{\mathsf{d}} \end{bmatrix}_{n+1}^{k+1}=-\begin{bmatrix} \boldsymbol{\mathsf{r}}_{u}^{\mathrm{s}} \\ \\ \boldsymbol{\mathsf{r}}_{p}^{\mathrm{s}} \\ \\ \boldsymbol{\mathsf{r}}_{v}^{\mathrm{f}} \\ \\ \boldsymbol{\mathsf{r}}_{p}^{\mathrm{f}} \\ \\ \boldsymbol{\mathsf{r}}_{\lambda} \\ \\ \boldsymbol{\mathsf{r}}_{\mathit{\Lambda}} \\ \\ \boldsymbol{\mathsf{r}}_{d}\end{bmatrix}_{n+1}^{k} \end{aligned}\end{split}\]

Fluid-reduced-Solid Interaction (FrSI)

– Boundary term:

\[\begin{aligned} r_{\tilde{\text{s}}} := \int\limits_{\mathit{\Gamma}_{0}^{\text{f}\text{-}\tilde{\text{s}}}} h_{0} \left[\rho_{0,s}\frac{\partial\boldsymbol{v}}{\partial t}\cdot\delta\boldsymbol{v} + \tilde{\boldsymbol{P}}(\boldsymbol{u}_{\text{f}}(\boldsymbol{v}) + \tilde{\boldsymbol{u}}_{\mathrm{pre}},\boldsymbol{v}) : \tilde{\nabla}_{0}\delta\boldsymbol{v}\right]\mathrm{d}A \end{aligned}\]

– Discrete nonlinear system to solve in each time step \(n\):

(61)\[\begin{split}\begin{aligned} \boldsymbol{\mathsf{r}}_{n+1} = \begin{bmatrix} \boldsymbol{\mathsf{V}}_{v}^{\mathit{\Gamma}^\mathrm{T}}\boldsymbol{\mathsf{r}}_{v}(\boldsymbol{\mathsf{V}}_{v}^{\mathit{\Gamma}}\tilde{\boldsymbol{\mathsf{v}}},\boldsymbol{\mathsf{p}},\boldsymbol{\mathsf{\Lambda}},\boldsymbol{\mathsf{d}}) \\ \boldsymbol{\mathsf{r}}_{p}(\boldsymbol{\mathsf{p}},\boldsymbol{\mathsf{V}}_{v}^{\mathit{\Gamma}}\tilde{\boldsymbol{\mathsf{v}}},\boldsymbol{\mathsf{d}}) \\ \boldsymbol{\mathsf{r}}_{\mathit{\Lambda}}(\boldsymbol{\mathsf{\Lambda}},\boldsymbol{\mathsf{V}}_{v}^{\mathit{\Gamma}}\tilde{\boldsymbol{\mathsf{v}}},\boldsymbol{\mathsf{d}}) \\ \boldsymbol{\mathsf{r}}_{d}(\boldsymbol{\mathsf{d}},\boldsymbol{\mathsf{V}}_{v}^{\mathit{\Gamma}}\tilde{\boldsymbol{\mathsf{v}}}) \end{bmatrix}_{n+1} = \boldsymbol{\mathsf{0}} \end{aligned}\end{split}\]

– Discrete linear system to solve in each Newton iteration \(k\):

(62)\[\begin{split}\begin{aligned} \begin{bmatrix} \boldsymbol{\mathsf{V}}_{v}^{\mathit{\Gamma}^\mathrm{T}}\boldsymbol{\mathsf{K}}_{vv}\boldsymbol{\mathsf{V}}_{v}^{\mathit{\Gamma}} & \boldsymbol{\mathsf{V}}_{v}^{\mathit{\Gamma}^\mathrm{T}}\boldsymbol{\mathsf{K}}_{vp} & \boldsymbol{\mathsf{V}}_{v}^{\mathit{\Gamma}^\mathrm{T}}\boldsymbol{\mathsf{K}}_{v\mathit{\Lambda}} & \boldsymbol{\mathsf{V}}_{v}^{\mathit{\Gamma}^\mathrm{T}}\boldsymbol{\mathsf{K}}_{vd} \\ \\ \boldsymbol{\mathsf{K}}_{pv}\boldsymbol{\mathsf{V}}_{v}^{\mathit{\Gamma}} & \boldsymbol{\mathsf{K}}_{pp} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{pd} \\ \\ \boldsymbol{\mathsf{K}}_{\mathit{\Lambda}v}\boldsymbol{\mathsf{V}}_{v}^{\mathit{\Gamma}} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{\mathit{\Lambda}\mathit{\Lambda}} & \boldsymbol{\mathsf{K}}_{\mathit{\Lambda}d} \\ \\ \boldsymbol{\mathsf{K}}_{dv}\boldsymbol{\mathsf{V}}_{v}^{\mathit{\Gamma}} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{dd} \end{bmatrix}_{n+1}^{k}\begin{bmatrix} \Delta\tilde{\boldsymbol{\mathsf{v}}} \\ \\ \Delta\boldsymbol{\mathsf{p}} \\ \\ \Delta\boldsymbol{\mathsf{\Lambda}}\\ \\ \Delta\boldsymbol{\mathsf{d}} \end{bmatrix}_{n+1}^{k+1}=-\begin{bmatrix} \boldsymbol{\mathsf{V}}_{v}^{\mathit{\Gamma}^\mathrm{T}}\boldsymbol{\mathsf{r}}_{v} \\ \\ \boldsymbol{\mathsf{r}}_{p} \\ \\ \boldsymbol{\mathsf{r}}_{\mathit{\Lambda}} \\ \\ \boldsymbol{\mathsf{r}}_{d}\end{bmatrix}_{n+1}^{k} \end{aligned}\end{split}\]

Multiphase Fluid

– Primary variables: velocity \(\boldsymbol{v}\), pressure \(p\), phase \(\phi\), and (chem.) potential \(\mu\)

– conservative form of Navier-Stokes, with capillary force term:

(63)\[\begin{split}\begin{aligned} \frac{\partial(\rho\boldsymbol{v})}{\partial t} + \nabla\cdot(\rho(\boldsymbol{v}\otimes\boldsymbol{v})) &= \nabla \cdot \boldsymbol{\sigma} -\phi\nabla\mu + \hat{\boldsymbol{b}} &&\text{in} \; \mathit{\mathit{\Omega}}_t \times [0, T], \\ \frac{\partial\rho}{\partial t} + \nabla\cdot(\rho\boldsymbol{v}) &= 0 &&\text{in} \; \mathit{\mathit{\Omega}}_t \times [0, T],\\ \boldsymbol{v} &= \hat{\boldsymbol{v}} &&\text{on} \; \mathit{\mathit{\Gamma}}_t^{\mathrm{D}} \times [0, T],\\ \boldsymbol{t} = \boldsymbol{\sigma}\boldsymbol{n} &= \hat{\boldsymbol{t}} &&\text{on} \; \mathit{\mathit{\Gamma}}_t^{\mathrm{N}} \times [0, T],\\ \boldsymbol{v}(\boldsymbol{x},0) &= \hat{\boldsymbol{v}}_{0}(\boldsymbol{x}) &&\text{in} \; \mathit{\mathit{\Omega}}_t, \end{aligned}\end{split}\]

with the constitutive equation for the Cauchy stress tensor given by

\[\begin{aligned} \boldsymbol{\sigma} = -p \boldsymbol{I} + \eta \left(\nabla \boldsymbol{v} + (\nabla \boldsymbol{v})^{\mathrm{T}}\right) + \left(\zeta - \frac{2}{d}\eta\right)\nabla\cdot\boldsymbol{v} \end{aligned}\]

– conservative form of the Cahn-Hilliard equations:

(64)\[\begin{split}\begin{aligned} \frac{\partial \phi}{\partial t} + \nabla\cdot(\phi\boldsymbol{v}) &= - \nabla\cdot\boldsymbol{J} &&\text{in} \; \mathit{\mathit{\Omega}}_t \times [0, T], \\ \mu - \frac{\mathrm{d}\psi}{\mathrm{d}\phi} + \nabla\cdot(\kappa \nabla \phi) &= 0 &&\text{in} \; \mathit{\mathit{\Omega}}_t \times [0, T],\\ \phi(\boldsymbol{x},0) &= \hat{\phi}_{0}(\boldsymbol{x}) &&\text{in} \; \mathit{\mathit{\Omega}}_t, \end{aligned}\end{split}\]

with the constitutive equation for the diffusive flux

\[\begin{aligned} \boldsymbol{J} = -M(\phi)\nabla(\mu + \alpha p) \end{aligned}\]

where \(\alpha=\frac{\rho_1-\rho_2}{\rho_1+\rho_2}\)

We assume \(\phi \in [a, b]\) and define the normalized phase field variable

\[\begin{aligned} \chi(\phi) = \frac{\phi-a}{b-a} \end{aligned}\]

and fluid density as well as dynamic and bulk viscosity read

\[\begin{split}\begin{aligned} \rho(\chi(\phi)) &= \rho_1 (1 - \chi) + \rho_2 \chi, \\ \eta(\chi(\phi)) &= \eta_1 (1 - \chi) + \eta_2 \chi, \\ \zeta(\chi(\phi)) &= \zeta_1 (1 - \chi) + \zeta_2 \chi \end{aligned}\end{split}\]

Chemical potential:

\[\begin{aligned} \psi(\phi) = D_0 (a-\phi)^2 (b-\phi)^2 \end{aligned}\]

with \(D_0\) a bulk free-energy parameter

Mobility:

\[ \begin{align}\begin{aligned}\begin{aligned} M(\phi) = M_0 |(a-\phi)^{\gamma} (b-\phi)^{\gamma}| :label: mobility\\\end{aligned}\end{aligned}\end{align} \]

Variational form:

Find \((\boldsymbol{v},p,\phi,\mu) \in \boldsymbol{\mathcal{V}}_{v}^{h,D} \times \mathcal{V}_{p}^{h,D} \times \mathcal{V}_{\phi}^{h,D} \times \mathcal{V}_{\mu}^{h,D}\) such that

\[\begin{split}\begin{aligned} &\int\limits_{\mathit{\Omega}} \left(\frac{\partial(\rho\boldsymbol{v})}{\partial t} + \nabla\cdot(\rho(\boldsymbol{v}\otimes\boldsymbol{v}))\right) \cdot \delta\boldsymbol{v} \,\mathrm{d}v + \int\limits_{\Omega}\boldsymbol{\sigma} : \nabla\delta\boldsymbol{v}\,\mathrm{d}v + \int\limits_{\Omega}\phi \nabla\mu\cdot\delta\boldsymbol{v}\,\mathrm{d}v - \delta\mathcal{P}_{\mathrm{ext}}(\boldsymbol{v};\delta\boldsymbol{v}) = 0, \\ &\int\limits_{\mathit{\Omega}}\left(\frac{\partial\rho}{\partial t} + \nabla\cdot(\rho\boldsymbol{v})\right)\delta p \,\mathrm{d}v = 0, \\ &\int\limits_{\mathit{\Omega}} \left(\frac{\partial \phi}{\partial t} + \nabla\cdot(\phi\boldsymbol{v})\right) \delta \phi \, \mathrm{d}v - \int\limits_{\Omega} \boldsymbol{J} \cdot \nabla \delta \phi \, \mathrm{d}v \quad + \text{(Neumann BCs)} = 0, \\ &\int\limits_{\mathit{\Omega}} \mu \,\delta\mu \, \mathrm{d}v - \int\limits_{\Omega} \frac{\mathrm{d}\psi}{\mathrm{d}\phi} \delta\mu \,\mathrm{d}v - \int\limits_{\Omega} \kappa \nabla \phi \cdot \nabla \delta\mu \, \mathrm{d}v \quad + \text{(Neumann BCs)} = 0, \end{aligned}\end{split}\]

for all \((\delta\boldsymbol{v},\delta p,\delta\phi,\delta\mu) \in \boldsymbol{\mathcal{V}}_{v}^{h} \times \mathcal{V}_{p}^{h} \times \mathcal{V}_{\phi}^{h} \times \mathcal{V}_{\mu}^{h}\), with

\[\begin{aligned} \frac{\partial\rho}{\partial t} = \rho^{\prime}(\phi)\frac{\partial\phi}{\partial t} \end{aligned}\]

– Discrete nonlinear system to solve in each time step \(n\):

(65)\[\begin{split}\begin{aligned} \boldsymbol{\mathsf{r}}_{n+1} = \begin{bmatrix} \boldsymbol{\mathsf{r}}_{v}(\boldsymbol{\mathsf{v}},\boldsymbol{\mathsf{p}},\boldsymbol{\upphi},\boldsymbol{\upmu}) \\ \boldsymbol{\mathsf{r}}_{p}(\boldsymbol{\mathsf{p}},\boldsymbol{\mathsf{v}},\boldsymbol{\upphi}) \\ \boldsymbol{\mathsf{r}}_{\phi}(\boldsymbol{\upphi},\boldsymbol{\mathsf{v}},\boldsymbol{\mathsf{p}},\boldsymbol{\upmu}) \\ \boldsymbol{\mathsf{r}}_{\mu}(\boldsymbol{\upmu},\boldsymbol{\upphi}) \end{bmatrix}_{n+1} = \boldsymbol{\mathsf{0}} \end{aligned}\end{split}\]

– Discrete linear system to solve in each Newton iteration \(k\):

(66)\[\begin{split}\begin{aligned} \begin{bmatrix} \boldsymbol{\mathsf{K}}_{vv} & \boldsymbol{\mathsf{K}}_{vp} & \boldsymbol{\mathsf{K}}_{v\phi} & \boldsymbol{\mathsf{K}}_{v\mu} \\ \\ \boldsymbol{\mathsf{K}}_{pv} & \boldsymbol{\mathsf{K}}_{pp} & \boldsymbol{\mathsf{K}}_{p\phi} & \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}\\ \\ \boldsymbol{\mathsf{K}}_{\phi v} & \boldsymbol{\mathsf{K}}_{\phi p} & \boldsymbol{\mathsf{K}}_{\phi \phi} & \boldsymbol{\mathsf{K}}_{\phi \mu} \\ \\ \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \textcolor{lightgray}{\boldsymbol{\mathsf{0}}}& \boldsymbol{\mathsf{K}}_{\mu\phi} & \boldsymbol{\mathsf{K}}_{\mu\mu} \end{bmatrix}_{n+1}^{k}\begin{bmatrix} \Delta\boldsymbol{\mathsf{v}} \\ \\ \Delta\boldsymbol{\mathsf{p}} \\ \\ \Delta\boldsymbol{\upphi} \\ \\ \Delta\boldsymbol{\upmu}\end{bmatrix}_{n+1}^{k+1}=-\begin{bmatrix} \boldsymbol{\mathsf{r}}_{v} \\ \\ \boldsymbol{\mathsf{r}}_{p} \\ \\ \boldsymbol{\mathsf{r}}_{\phi} \\ \\ \boldsymbol{\mathsf{r}}_{\mu}\end{bmatrix}_{n+1}^{k} \end{aligned}\end{split}\]

Demos

Demo: Solid

– Physics description given in sec. 4.1

– Input files: demos/solid/

Cantilever under tip load

This example demonstrates how to set up a quasi-static solid mechanics elasticity problem. The deformation of a steel cantilever under transverse conservative load is simulated. The structure is fixed on one end. Quadratic 27-node hexahedral finite elements are used for the discretization of the domain. The well-known St. Venant-Kirchhoff material is used as constitutive law, which is a generalization of Hooke’s law to the nonlinear realm.

Study the setup shown in fig. 1 and the comments in the input file solid_cantilever.py Run the simulation, either in one of the provided Docker containers or using your own FEniCSx/Ambit installation, using the command

mpiexec -n 1 python3 solid_cantilever.py

It is fully sufficient to use one core (mpiexec -n 1) for the presented setup.

Open the results file results_solid_cantilever_displacement.xdmf in Paraview, and visualize the deformation over time.

Figure 2 shows the displacement magnitude at the end of the simulation.

Demo: Fluid

– Physics description given in sec. 4.2

– Input files: demos/fluid

2D channel flow

This example shows how to set up 2D fluid flow in a channel around a rigid obstacle. Incompressible Navier-Stokes flow is solved using Taylor-Hood elements (9-node biquadratic quadrilaterals for the velocity, 4-node bilinear quadrilaterals for the pressure).

Study the setup and the comments in the input file fluid_channel.py. Run the simulation, either in one of the provided Docker containers or using your own FEniCSx/Ambit installation, using the command

mpiexec -n 1 python3 fluid_channel.py

It is fully sufficient to use one core (mpiexec -n 1) for the presented setup.

Open the results file results_fluid_channel_velocity.xdmf and

results_fluid_channel_pressure.xdmf in Paraview, and visualize the velocity as well as the pressure over time.

Fig. 4 shows the velocity magnitude (top) as well as the pressure (bottom part) at the end of the simulation.

Demo: 0D flow

– Physics description given in sec. 4.3

– Input files: demos/flow0d

Systemic and pulmonary circulation

This example demonstrates how to simulate a cardiac cycle using a lumped-parameter (0D) model for the heart chambers and the entire circulation. Multiple heart beats are run until a periodic state criterion is met (which compares variable values at the beginning to those at the end of a cycle, and stops if the relative change is less than a specified value, here `eps_periodic' in the TIME_PARAMS dictionary). The problem is set up such that periodicity is reached after 5 heart cycles.

Study the setup in fig. 5 and the comments in the input file flow0d_heart_cycle.py. Run the simulation, either in one of the provided Docker containers or using your own FEniCSx/Ambit installation, using the command

python3 flow0d_heart_cycle.py

For postprocessing of the time courses of pressures, volumes, and fluxes of the 0D model, either use your own tools to plot the text output files (first column is time, second is the respective quantity), or make sure to have Gnuplot (and TeX) installed and navigate to the output folder (tmp/) in order to execute the script flow0d_plot.py (which lies in ambit/src/ambit_fe/postprocess/):

flow0d_plot.py -s flow0d_heart_cycle -n 100

A folder plot_flow0d_heart_cycle is created inside tmp/. Look at the results of pressures (\(p\)), volumes (\(V\)), and fluxes (\(q\), \(Q\)) over time. Subscripts v, at, ar, ven refer to ‘ventricular’, ‘atrial’, ‘arterial’, and ‘venous’, respectively. Superscripts l, r, sys, pul refer to ‘left’, ‘right’, ‘systemic’, and ‘pulmonary’, respectively. Try to understand the time courses of the respective pressures, as well as the plots of ventricular pressure over volume. Check that the overall system volume is constant and around 4-5 liters.

The solution is depicted in fig. 6, showing the time course of volumes and pressures of the circulatory system.

Demo: Solid + 0D flow

– Physics description given in sec. 4.4.1

– Input files: demos/solid_flow0d

3D heart, coupled to systemic and pulmonary circulation

This example demonstrates how to set up and simulate a two-chamber (left and right ventricular) solid mechanics heart model coupled to a closed-loop 0D circulatory system. A full dynamic heart cycle of duration 1 s is simulated, where the active contraction is modeled by a prescribed active stress approach. Passive material behavior of the heart muscle is governed by the Holzapfel-Ogden anisotropic strain energy function [HO09] and a strain rate-dependent viscous model [CTMS12]. We start the simulation with “prestressing” using the MULF method [GForsterW10, SG21], which allows to imprint loads without changing the geometry, where the solid is loaded to the initial left and right ventricular pressures. Thereafter, we kickstart the dynamic simulation with passive ventricular filling by the systole of the atria (0D chamber models). Ventricular systole happens in \(t \in [0.2\;\mathrm{s}, 0.53\;\mathrm{s}]\), hence lasting a third of the whole cycle time. After systole, the heart relaxes and eventually fills to about the same pressure as it has been initialized to.

NOTE: For demonstrative purposes, a fairly coarse finite element discretization is chosen here, which by no means yields a spatially converged solution and which may be prone to locking phenomena. The user may increse the parameter `order_disp' in the FEM_PARAMS section from 1 to 2 (and increase `quad_degree' to 6) such that quadratic finite element ansatz functions (instead of linear ones) are used. While this will increase accuracy and mitigate locking, computation time will increase.

Study the setup shown in fig. 7 and the comments in the input file solid_flow0d_heart_cycle.py. Run the simulation, either in one of the provided Docker containers or using your own FEniCSx/Ambit installation, using the command

mpiexec -n 1 python3 solid_flow0d_heart_cycle.py

It is fully sufficient to use one core (mpiexec -n 1) for the presented setup, while you might want to use more (e.g., mpiexec -n 4) if you increase `order_disp' to 2.

Open the results file results_solid_flow0d_heart_cycle_displacement.xdmf in Paraview, and visualize the deformation over the heart cycle.

For postprocessing of the time courses of pressures, volumes, and fluxes of the 0D model, either use your own tools to plot the text output files (first column is time, second is the respective quantity), or make sure to have Gnuplot (and TeX) installed and navigate to the output folder (tmp/) in order to execute the script flow0d_plot.py (which lies in ambit/src/ambit_fe/postprocess/):

flow0d_plot.py -s solid_flow0d_heart_cycle -V0 117e3 93e3 0 0 0

A folder plot_solid_flow0d_heart_cycle is created inside tmp/. Look at the results of pressures (\(p\)), volumes (\(V\)), and fluxes (\(q\), \(Q\)) over time. Subscripts v, at, ar, ven refer to ‘ventricular’, ‘atrial’, ‘arterial’, and ‘venous’, respectively. Superscripts l, r, sys, pul refer to ‘left’, ‘right’, ‘systemic’, and ‘pulmonary’, respectively. Try to understand the time courses of the respective pressures, as well as the plots of ventricular pressure over volume. Check that the overall system volume is constant and around 4-5 liters.

NOTE: This setup computes only one cardiac cycle, which does not yield a periodic state solution (compare e.g. initial and end-cyclic right ventricular pressures and volumes, which do not coincide). Change the parameter number_of_cycles from 1 to 10 and re-run the simulation. The simulation will stop when the cycle error (relative change in 0D variable quantities from beginning to end of a cycle) falls below the value of `eps_periodic' (set to \(5 \%\)). How many cycles are needed to reach periodicity?

Figure 8 shows a high-fidelity solution using a refined mesh and quadratic tetrahedral elements. Compare your solution from the coarser mesh. What is the deviation in ventricular volume?

Demo: Fluid + 0D flow

– Physics description given in sec. 4.4.2

– Input files: demos/fluid_flow0d

Blocked pipe flow with 0D model bypass

This example demonstrates how to couple 3D fluid flow to a 0D lumped-parameter model. Incompressible transient Navier-Stokes flow in a pipe with prescribed inflow is solved, with the special constraint that an internal boundary (all-time closed valve) separates region 1 and region 2 of the pipe. This internal Dirichlet condition can only be achieved by splitting the pressure space, hence having duplicate pressure nodes at the valve plane. Otherwise, fluid would experience deceleration towards the valve and unphysical acceleration behind it, since the pressure gradient drives fluid flow.

This example demonstrates how the closed valve can be bypassed by a 0D flow model that links the 3D fluid out-flow of one region to the in-flow of the other region. The 0D model consists of two Windkessel models in series, each having compliance, resistance, and inertance elements.

Study the setup shown in fig. 9 and the comments in the input file fluid_flow0d_pipe.py. Run the simulation, either in one of the provided Docker containers or using your own FEniCSx/Ambit installation, using the command

mpiexec -n 1 python3 fluid_flow0d_pipe.py

It is fully sufficient to use one core (mpiexec -n 1) for the presented setup.

Open the results file results_fluid_flow0d_pipe_velocity.xdmf in Paraview, and visualize the velocity over time.

Think of which parameter(s) of the 0D model to tweak in order to achieve a) little to no fluid in-flow (into \(\mathit{\Gamma}_{\mathrm{in}}^{\mathrm{f\text{-}0d}}\)), b) almost the same flow across \(\mathit{\Gamma}_{\mathrm{out}}^{\mathrm{f\text{-}0d}}\) and \(\mathit{\Gamma}_{\mathrm{in}}^{\mathrm{f\text{-}0d}}\). Think of where the flow is going to in case of a).

Figure shows the velocity streamlines and magnitude at the end of the simulation.

Demo: FSI

– Physics description given in sec. 4.4.4

– Input files: demos/fsi

Channel flow around elastic flag

Incompressible fluid flow in a 2D channel around an elastic flag is studied. The setup corresponds to the well-known Turek benchmark [TH06]. Here, the two cases FSI2 and FSI3 from the original setup are investigated. A prescribed inflow velocity with parabolic inflow profile is used:

\[\begin{aligned} \boldsymbol{v}_{\text{f}}= \bar{v}(t,y) \boldsymbol{e}_{x} \quad & \text{on}\; \mathit{\Gamma}_{t,\mathrm{in}}^{D,\text{f}} \times [0,T], \label{eq:flag_dbc_in} \end{aligned}\]

with

\[\begin{split}\begin{aligned} \bar{v}(t,y) = \begin{cases} 1.5 \,\bar{U}\, \frac{y(H-y)}{\left(\frac{H}{2}\right)^2} \frac{1-\cos\left(\frac{\pi}{2}t\right)}{2}, & \text{if} \; t < 2, \\ 1.5 \,\bar{U}\, \frac{y(H-y)}{\left(\frac{H}{2}\right)^2}, & \text{else}, \end{cases} \label{eq:flag_dbcs_func} \end{aligned}\end{split}\]

with \(\bar{U}=10^{3}\;\mathrm{mm}/\mathrm{s}\) (FSI2) and \(\bar{U}=2\cdot 10^{3}\;\mathrm{mm}/\mathrm{s}\) (FSI3).

Geometrical parameters, given in \([\mathrm{mm}]\), are:

\(L\)	\(H\)	\(r\)	\(l\)	\(h\)	\(d_x\)	\(d_y\)
2500	410	50	350	20	200	200

The fluid is discretized with quadrilateral \(\mathbb{Q}^2\)-\(\mathbb{Q}^1\) Taylor-Hood finite elements, hence no stabilization is needed. The solid is discretized with displacement-based quadrilateral \(\mathbb{Q}^2\) elements. Temporal discretization for both the solid and the fluid are carried out with a Generalized-\(\alpha\) scheme with no numerical damping (\(\rho_{\mathrm{inf}}=1\)).

Study the setup shown in fig. 11 together with the parameters in the table and the comments in the input file fsi_channel_flag.py. Run the simulation for FSI2 and FSI3 cases, either in one of the provided Docker containers or using your own FEniCSx/Ambit installation, using the command

mpiexec -n 1 python3 fsi_channel_flag.py

If your system allows, use more than one core (e.g. mpiexec -n 4) in order to speed up the simulation a bit.

The physics of the problem are strongly time-dependent, and a (near-)periodic oscillation of the flag only occurs after \(t\approx 10\;\mathrm{s}\) (FSI2) and \(t\approx 5\;\mathrm{s}\) (FSI3). Run the problem to the end (\(t = 15\;\mathrm{s}\) for FSI2, \(t = 7.5\;\mathrm{s}\) for FSI3), be patient, and monitor the flag tip displacement over time.

Figure 12 depicts the velocity at three instances in time towards the end of the simulation for the FSI2 case, and figure 13 shows the flag’s tip displacement over time compared to the reference solution, over a time interval where the solution has become periodic. (Note that the reference solution from the official link shown in the input file needs to be time-adjusted, i.e. synchronized with the first peak in the interval of interest, since the time column of the reference data does not correspond to the physical time of the problem setup.)

Demo: Multiphase fluid

– Physics description given in sec. 4.4.7

– Input files: demos/fluid_multiphase

Rising bubble under gravity

We study a 2D multiphase flow benchmark, where two incompressible fluids of different densities and different viscosities interact: a fluid of lower density (and viscosity) is immersed into a surrounding fluid of higher density (and viscosity), and a gravitational force acts in vertical direction, cf. the setup in Fig. 14. This setup corresponds to the benchmark case described in [BtE26],:cite:p:ten-eikelder2024. The following Dirichlet boundary conditions apply:

\[\begin{split}\begin{aligned} \boldsymbol{v} = \boldsymbol{0} \quad & \text{on}\; \mathit{\Gamma}_{\mathrm{b}}^{D}\cup\mathit{\Gamma}_{\mathrm{t}}^{D} \times [0,T], \\ \boldsymbol{v}\cdot\boldsymbol{n} = 0 \quad & \text{on}\; \mathit{\Gamma}_{\mathrm{r}}^{D}\cup\mathit{\Gamma}_{\mathrm{l}}^{D} \times [0,T]. \label{eq:rising_dbc} \end{aligned}\end{split}\]

The body force term takes the form \(\hat{\boldsymbol{b}}=-\rho\boldsymbol{g}\), with \(\boldsymbol{g}\cdot\boldsymbol{e}_{y}=0.98\). Our phase field variable is bounded between \(a=-1\) and \(b=1\). Further, Cahn-Hilliard parameters are chosen as

\[\begin{aligned} D_0 = \frac{\tilde{\sigma}}{4\epsilon}, \qquad M_0 = 0.1\epsilon^2, \qquad \kappa = \tilde{\sigma}\epsilon, \end{aligned}\]

with \(\tilde{\sigma}=\frac{3\sigma}{2\sqrt{2}}\) (surface energy density coefficient) and \(\epsilon=0.64\,h_{\mathrm{e}}\), where \(h_{\mathrm{e}}\) is the finite element edge length. The mobility exponent in Eq. ([equation-mobility]) is chosen to \(\gamma=1\). Parameters for the two cases are shown in the following table. Bulk viscosities are assumed zero (\(\zeta=0\)).

	math: rho_1	math: rho_2	math: eta_1	math: eta_2	math: sigma
Case 1	1000	100	10	1	24.5
Case 2	1000	1	10	0.1	1.96

The fluid is discretized with quadrilateral \(\mathbb{Q}^2\)-\(\mathbb{Q}^1\) Taylor-Hood finite elements, while a quadrilateral \(\mathbb{Q}^1\)-\(\mathbb{Q}^1\) discretization is chosen for the Cahn-Hilliard equations. A BDF2 time-integration scheme is chosen for both fluid and Cahn-Hilliard.

Study the setup shown in Fig. 14 together with the parameters in the table and the comments in the input file fluid_multiphase_rising_bubble.py. Run the simulation for both cases 1 and 2, either in one of the provided Docker containers or using your own FEniCSx/Ambit installation, using the command

mpiexec -n 1 python3 fluid_multiphase_rising_bubble.py

Try different settings for the mesh, e.g. vary meshsize inside mesh_domain (increase from [32,64] to [64,128] to [128,256]; the first value is the number of elements used in horizontal, the second the number in vertical direction). If your system allows, use more than one core, especially for finer mesh sizes (hence run with mpiexec -n 4 ... or more). A good rule of thumb is to have 10000–20000 degrees of freedom per process. The time step and the interface thickness parameter \(\epsilon\) will be changed automatically depending on the chosen mesh size. Compare the results for the two cases and the different mesh sizes.

Figure 15 shows the results of phase field and magnitude of fluid velocity for a mesh size of \(128 \times 256\) for both cases at the end of the simulation time \(t=3\).

Figure 16 plots the center of mass as well as the rise velocity of the lower-density fluid bubble over time (mesh size \(128 \times 256\)) and compares the solution with two different reference solutions (one a mass-, the other a volume-averaged velocity formulation) from Brunk and ten Eikelder [BtE26] as well as ten Eikelder and Schillinger [tES24]. We observe good agreements between our (mass-averaged velocity) solution and the volume-averaged velocity reference solutions. Slight deviations to the mass-averaged reference solutions are observed.

Table of symbols

\[\begin{split}\nonumber \begin{aligned} &\mathit{\Omega}_0,\mathit{\Omega}&&: \text{reference, current domain} \\ &\mathit{\Gamma}_0,\mathit{\Gamma}&&: \text{reference, current boundary} \\ &\boldsymbol{x}_0, \boldsymbol{x} &&: \text{coordinates of the reference, current frame} \\ &\boldsymbol{e}_x, \boldsymbol{e}_y, \boldsymbol{e}_z &&: \text{unit vectors of the cartesian reference frame} \\ &\boldsymbol{n}_0, \boldsymbol{n} &&: \text{unit outward normal defined in the reference, current frame} \\ &\nabla_{0},\nabla &&: \text{Nabla operator with respect to the reference, current frame} \\ &\nabla\boldsymbol{v} := \frac{\partial v_i}{\partial x_j} \boldsymbol{e}_{i} \otimes\boldsymbol{e}_{j} &&: \text{Gradient of a vector field} \\ &\nabla\cdot\boldsymbol{v} := \frac{\partial v_i}{\partial x_j} \boldsymbol{e}_{i} \cdot\boldsymbol{e}_{j} &&: \text{Divergence of a vector field} \\ &\boldsymbol{e}_{i} &&: \text{Basis vectors}, i\in \{1,2,3\} \\ &t, T &&: \text{current, end time of an initial (boundary) value problem} \\ &\boldsymbol{u}, \hat{\boldsymbol{u}}_{0} &&: \text{solid mechanics displacement field, and prescribed initial value} \\ &\delta\boldsymbol{u}, \Delta\boldsymbol{u} &&: \text{solid mechanics displacement test, trial function} \\ & p &&: \text{solid mechanics hydrostatic pressure, or fluid mechanics pressure} \\ & \delta p, \Delta p &&: \text{solid or fluid mechanics pressure test, trial function} \\ &\boldsymbol{v}=\frac{\partial\boldsymbol{u}}{\partial t}, \hat{\boldsymbol{v}}_{0} &&: \text{solid mechanics velocity, and prescribed initial value} \\ &\boldsymbol{a}=\frac{\partial^2\boldsymbol{u}}{\partial t^2} &&: \text{solid mechanics acceleration} \\ &\boldsymbol{v}, \hat{\boldsymbol{v}}_{0} &&: \text{fluid mechanics velocity, and prescribed initial value} \\ &\delta\boldsymbol{v}, \Delta\boldsymbol{v} &&: \text{fluid mechanics velocity test, trial function} \\ &\boldsymbol{a}=\frac{\partial\boldsymbol{v}}{\partial t} &&: \text{fluid mechanics acceleration} \\ &\boldsymbol{d} &&: \text{ALE domain displacement} \\ &\delta\boldsymbol{d}, \Delta\boldsymbol{d} &&: \text{ALE domain displacement test, trial function} \\ &\hat{\boldsymbol{b}}_0, \hat{\boldsymbol{b}} &&: \text{body force vector defined in the reference, current frame} \\ &\widehat{\boldsymbol{w}}=\frac{\partial\boldsymbol{d}}{\partial t} &&: \text{ALE domain velocity} \\ &\rho_0, \rho &&: \text{reference, current density} \\ &\phi, \mu &&: \text{phase field, potential in Cahn-Hilliard equations} \end{aligned}\end{split}\]

\[\begin{split}\nonumber \begin{aligned} &\boldsymbol{P}=\boldsymbol{F}\boldsymbol{S} &&: \text{1st Piola Kirchhoff stress tensor} \\ &\boldsymbol{F}=\boldsymbol{I}+\nabla_{0}\boldsymbol{u} &&: \text{solid deformation gradient} \\ &\widehat{\boldsymbol{F}}=\boldsymbol{I}+\nabla_{0}\boldsymbol{d} &&: \text{ALE deformation gradient} \\ &J=\det \boldsymbol{F} &&: \text{determinant of solid deformation gradient} \\ &\widehat{J}=\det \widehat{\boldsymbol{F}} &&: \text{determinant of ALE deformation gradient} \\ &\boldsymbol{S} &&: \text{2nd Piola-Kirchhoff stress tensor} \\ &\boldsymbol{\sigma} &&: \text{Cauchy stress tensor} \\ &\boldsymbol{t}_0, \hat{\boldsymbol{t}}_{0} &&: \text{1st Piola-Kirchhoff traction, prescribed 1st Piola-Kirchhoff traction} \\ &\boldsymbol{t}, \hat{\boldsymbol{t}} &&: \text{Cauchy traction, prescribed Cauchy traction} \end{aligned}\end{split}\]

[ALA+16]

C. J. Arthurs, K. D. Lau, K. N. Asrress, S. R. Redwood, and C. A. Figueroa. A mathematical model of coronary blood flow control: simulation of patient-specific three-dimensional hemodynamics during exercise. Am J Physiol Heart Circ Physiol, 310(9):H1242–H1258, 2016. doi:10.1152/ajpheart.00517.2015.

[BAA+22]

S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, and J. Zhang. PETSc/TAO users manual. Technical Report ANL-21/39 - Revision 3.17, Argonne National Laboratory, 2022.

[BtE26] (1,2,3)

A. Brunk and M. F. P. ten Eikelder. A simple, fully-discrete, unconditionally energy-stable method for the two-phase Navier-Stokes Cahn-Hilliard model with arbitrary density ratios. Journal of Computational Physics, 548():114558, 2026. doi:10.1016/j.jcp.2025.114558.

[CTMS12]

D. Chapelle, P. Le Tallec, P. Moireau, and M. Sorine. Energy-preserving muscle tissue model: formulation and compatible discretizations. Journal for Multiscale Computational Engineering, 10(2):189–211, 2012. doi:10.1615/IntJMultCompEng.2011002360.

[FACC14]

C. Farhat, P. Avery, T. Chapman, and J. Cortial. Dimensional reduction of nonlinear finite element dynamic models with finite rotations and energy-based mesh sampling and weighting for computational efficiency. International Journal for Numerical Methods in Engineering, 98(9):625–662, 2014. doi:10.1002/nme.4668.

[GForsterW10] (1,2)

M. W. Gee, Ch. Förster, and W. A. Wall. A computational strategy for prestressing patient-specific biomechanical problems under finite deformation. International Journal for Numerical Methods in Biomedical Engineering, 26(1):52–72, 2010. doi:10.1002/cnm.1236.

[GKL09]

M. W. Gee, C. T. Kelley, and R. B. Lehoucq. Pseudo-transient continuation for nonlinear transient elasticity. International Journal for Numerical Methods in Engineering, 78(10):1209–1219, 2009. doi:10.1002/nme.2527.

[GCM95]

J. M. Guccione, K. D. Costa, and A. D. McCulloch. Finite element stress analysis of left ventricular mechanics in the beating dog heart. J Biomech, 28(10):1167–1177, 1995. doi:10.1016/0021-9290(94)00174-3.

[GoktepeAPK10]

S. Göktepe, O. J. Abilez, K. K. Parker, and E. Kuhl. A multiscale model for eccentric and concentric cardiac growth through sarcomerogenesis. J Theor Biol, 265(3):433–442, 2010. doi:10.1016/j.jtbi.2010.04.023.

[Hir19] (1,2,3)

M. Hirschvogel. Computational modeling of patient-specific cardiac mechanics with model reduction-based parameter estimation and applications to novel heart assist technologies. Verlag Dr. Hut, MediaTUM, 1 edition, 2019. URL: https://mediatum.ub.tum.de/1445317.

[HBBN24]

M. Hirschvogel, M. Balmus, M. Bonini, and D. Nordsletten. Fluid-reduced-solid interaction (FrSI): physics- and projection-based model reduction for cardiovascular applications. Journal of Computational Physics, 506:112921, 2024. URL: https://www.sciencedirect.com/science/article/pii/S0021999124001700, doi:10.1016/j.jcp.2024.112921.

[HBJ+17]

M. Hirschvogel, M. Bassilious, L. Jagschies, S. M. Wildhirt, and M. W. Gee. A monolithic 3D-0D coupled closed-loop model of the heart and the vascular system: Experiment-based parameter estimation for patient-specific cardiac mechanics. Int J Numer Method Biomed Eng, 33(8):e2842, 2017. doi:10.1002/cnm.2842.

[Hir24]

Marc Hirschvogel. Ambit – A FEniCS-based cardiovascular multi-physics solver. Journal of Open Source Software, 9(93):5744, 2024. URL: https://doi.org/10.21105/joss.05744, doi:10.21105/joss.05744.

[Hol00]

G. A. Holzapfel. Nonlinear Solid Mechanics – A Continuum Approach for Engineering. Wiley Press Chichester, 2000.

[HO09] (1,2)

G. A. Holzapfel and R. W. Ogden. Constitutive modelling of passive myocardium: A structurally based framework for material characterization. Phil Trans R Soc A, 367(1902):3445–3475, 2009. doi:10.1098/rsta.2009.0091.

[LMW12]

A. Logg, K.-A. Mardal, and G. N. Wells, editors. Automated Solution of Differential Equations by the Finite Element Method – The FEniCS Book. Springer, 2012. ISBN 978-3-642-23098-1. doi:10.1007/978-3-642-23099-8.

[NMK+11]

D. A. Nordsletten, M. McCormick, P. J. Kilner, P. Hunter, D. Kay, and N. P. Smith. Fluid-solid coupling for the investigation of diastolic and systolic human left ventricular function. International Journal for Numerical Methods in Biomedical Engineering, 27(7):1017–1039, 2011. doi:10.1002/cnm.1405.

[SG21] (1,2)

A. Schein and M. W. Gee. Greedy maximin distance sampling based model order reduction of prestressed and parametrized abdominal aortic aneurysms. Advanced Modeling and Simulation in Engineering Sciences, 8(18):, 2021. doi:10.1186/s40323-021-00203-7.

[tES24]

M. F. P. ten Eikelder and D. Schillinger. The divergence-free velocity formulation of the consistent Navier-Stokes Cahn-Hilliard model with non-matching densities, divergence-conforming discretization, and benchmarks. Journal of Computational Physics, 513():113148, 2024. doi:10.1016/j.jcp.2024.113148.

[TO00] (1,2)

T. E. Tezduyar and Y. Osawa. Finite element stabilization parameters computed from element matrices and vectors. Computer Methods in Applied Mechanics and Engineering, 190(3–4):411–430, 2000. doi:10.1016/S0045-7825(00)00211-5.

[TH06]

S. Turek and J. Hron. Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow. In H.-J.Bungartz and M. Schäfer, editors, Lecture Notes in Computational Science and Engineering, volume 53, pages 371–385. Springer Berlin Heidelberg, 2006. doi:10.1007/3-540-34596-5\_15.

[WLW09]

N. Westerhof, J.-W. Lankhaar, and B. E. Westerhof. The arterial Windkessel. Med Biol Eng Comput, 47(2):H81–H88, 2009. doi:10.1007/s11517-008-0359-2.