svMultiPhysics - Beta Version

Introduction

svMultiPhysics is an open-source, parallel, finite element multi-physics solver providing capabilities to simulate the partial differential equations (PDEs) governing solid and fluid mechanics, diffusion, and electrophysiology. Equations can be solved coupled to simulate the interaction between multiple regions representing different physical systems. For example, in a coupled fluid-solid simulation the motion of the fluid can deform a solid region while the changing geometry of the solid region change the way the fluid flows. Equation coupling provides a framework for the computational modeling of whole heart dynamics.

svMultiPhysics implements boundary conditions useful for the simulation of blood flow in vascular models. It also has an interface to the svZeroDSolver that can be used for custom lumped parameter network (LPN) boundary conditions.

The svMultiPhysics solver is implemented in C++, a popular general-purpose object-oriented programming language that provides powerful features for producing modular software with a clear structure. C++ has many built-in data structures that can be used to efficiently organize and store data. The svMultiPhysics code is essentially a direct line-by-line translation of the Fortran svFSI solver code into C++.

The following sections describe the solver modeling capabilities, input data and format of the solver input parameters file used to run a simulation.

Documentation describing building, installing and testing svMultiPhysics can be found using the following links
Installing and building the solver
Testing Guide

Overview

The following sections provide an overview of the capabilities of the svMultiPhysics solver. In addition, they will provide a background for concepts and terminology that will be useful in later more detailed discussions of the solver.

Modeling capabilities

svMultiPhysics can numerically solve PDEs governing solid and fluid mechanics, diffusion, and electrophysiology. Some equations can be coupled and solved simultaneously to simulate the interaction between multiple regions representing different physical system.

svMultiPhysics supports the solution of the following equations

Boundary conditions

Boundary conditions are an integral part of a problem definition and act as constraints necessary for the solution of a PDE.

Different types of boundary conditions can be imposed on the boundary of the domain

• Dirichlet boundary condition - Fixed solution values
• Neumann boundary condition - Derivative of solution values
• Robin boundary condition - Weighted combination of Dirichlet and Neumann boundary conditions

A temporal and spatial distribution of values can be applied to the boundary of a domain.

svMultiPhysics supports boundary conditions useful for the simulation of blood flow in vascular models

• Windkessel (RCR) boundary condition - Electric circuit analogue that has a proximal resistance R in series with a parallel arrangement of a capacitance C and a distal resistance Rd
• Custom lumped parameter network (LPN) boundary condition - Implemented using an interface to the svZeroDSolver

Finite element method

The finite element method (FEM) is used to numerically solve transient PDEs governing the behavior of a physical system in two or three space dimensions. FEM approximates the geometry of a physical domain by subdividing it into a collection of discrete finite elements called a finite element mesh. The finite elements use numerical interpolation to approximate field variables (e.g. velocity) within a geometric region of a physical system. Elements can be implemented using a combination of linear, quadratic, and cubic interpolating polynomials.

svMultiPhysics supports the following element types

• line - linear and quadratic
• triangle - linear and quadratic
• quadrilateral - bilinear, serendipity, biquadratic
• tetrahedron - linear and quadratic
• hexahedron - trilinear, quadratic/serendipity, triquadratic
• wedge - linear

The finite element mesh together with appropriate physical properties and boundary conditions generate a system of algebraic equations. Depending on the PDE being solved this system of equations can be linear or nonlinear. The solution of a linear system can be directly solved to produce an approximation to the actual solution to the PDE. A nonlinear system must be solved using an iterative procedure to generate and solve a series of linear systems of equations converging to an approximation to the actual solution to the nonlinear PDE.

Most automatic mesh generation software packages produce finite element meshes are composed primarily of tetrahedron.

Numerical linear algebra

Numerical linear algebra uses computer algorithms to solve the linear system generated by FEM. Linear algebra libraries provide access to specialized or general purpose routines implementing a significant number of computer algorithms useful when solving linear systems.

svMultiPhysics supports interfaces to the following numerical linear algebra packages

• Custom numerical linear algebra routines included in the svMultiPhysics implementation
• Portable, Extensible Toolkit for Scientific Computation (PETSc)
• Trilinos

Parallel computation

The system of equations generated by FEM can often be quite large, involving millions of unknowns for the accurate representation of the vasculature of an anatomical region. The numerical solution of such a large system typically requires it to be distributed onto on a high-performance computational (HPC) cluster and solved in parallel. The svMultiPhysics solver uses Message Passing Interface (MPI) library routines to automatically partition and distribute an FEM simulation to run in parallel on an HPC cluster.

Mathematical preliminaries

The following use index (Einstein) notation.

Fluid mechanics

Strong form

The incompressible Navier-Stokes equations governing fluid flow are

$$ \rho\left(\frac{d\boldsymbol{u}}{dt} + \boldsymbol{u} \cdot \boldsymbol{\nabla} \boldsymbol{u} - \boldsymbol{b}\right) = \boldsymbol{\nabla} \cdot \boldsymbol{\sigma}, $$

$$ \boldsymbol{\nabla} \cdot \boldsymbol{u} = 0, $$

where $\boldsymbol{u} = \boldsymbol{u}\left(\boldsymbol{x}, t\right)$ is the velocity, $p = p\left(\boldsymbol{x}, t\right)$ is the pressure, $\boldsymbol{b} = \boldsymbol{b}\left(\boldsymbol{x}, t\right)$ is the body force, and $\rho$ is the fluid density. The first equation corresponds to the momentum conservation in the flow and the second equation corresponds to mass conservation. The momentum equation can augmented with a Darcy permeability term, $-\frac{\mu}{K}\boldsymbol{u}$, on the right-hand side to yield the Navier-Stokes-Brinkman equation [6] [7], such that

$$ \rho\left(\frac{d\boldsymbol{u}}{dt} + \boldsymbol{u} \cdot \boldsymbol{\nabla} \boldsymbol{u} - \boldsymbol{b}\right) = \boldsymbol{\nabla} \cdot \boldsymbol{\sigma} - \frac{\mu}{K}\boldsymbol{u}, $$

This equation models incompressible fluid flow in porous media. Here $K$ is the permeability of the porous media. With this equation, we can of course recover the Navier-Stokes equations by simply removing the Darcy component (i.e., $K \rightarrow \infty$).

The Cauchy stress tensor is $\boldsymbol{\sigma} = \boldsymbol{\sigma}\left(\boldsymbol{x}, t\right) = -p\boldsymbol{I} + 2\mu\left(\boldsymbol{u}\right)\epsilon$, where $\epsilon = \epsilon\left(\boldsymbol{u}\right) = \nabla^{s} \boldsymbol{u} = \frac{1}{2}\left(\nabla \boldsymbol{u} + \left(\nabla\boldsymbol{u}\right)^{\text{T}} \right)$ is the strain rate tensor. The effective dynamic viscosity, $\mu\left(\boldsymbol{u}\right)$, is written generally as a function of velocity here to account for non-Newtonian fluids. For Newtonian fluids, $\mu$ is a simply constant. The divergence of the Cauchy stress tensor, written in both vector and index notation, is

$$ \boldsymbol{\nabla} \cdot \boldsymbol{\sigma} = -\boldsymbol{\nabla}p + 2\boldsymbol{\epsilon}\boldsymbol{\nabla}\mu + \mu\nabla^{2}\boldsymbol{u}, $$

$$ \sigma_{ij,j} = -p_{,i} + 2\epsilon_{ij}\frac{\partial \mu}{\partial x_{j}} + \mu u_{i,kk}, $$

The boundary conditions are

$$ \boldsymbol{u} = \boldsymbol{g}, $$

$$ \boldsymbol{\sigma} \cdot \boldsymbol{n} = \boldsymbol{h}, $$

where $\boldsymbol{g}$ is the prescribed velocity and $\boldsymbol{h}$ is the prescribed traction.

We will solve the Navier-Stokes and Navier-Stokes-Brinkman equations numerically, using the finite element method for spatial discretization [8].

Standard (Galerkin) weak form

For the finite element method, we will first derive the Galerkin weak form for the Navier-Stokes-Brinkman equations. We define our trial and weighting function spaces,

$$ u_{i} \in \tau_{i} : \{u_{i} \in H^{1}\left(\Omega\right) \mid u_{i} = g_{i} \text{ on } \Gamma_{g_{i}}\}, $$

$$ w_{i} \in \nu_{i} : \{w_{i} \in H^{1}\left(\Omega\right) \mid w_{i} = 0 \text{ on } \Gamma_{g_{i}}\}, $$

$$ p, q \in Q : \{p \in L^{2}\left(\Omega\right)\}, $$

where $\boldsymbol{w}$ is the weighting function for velocity and $q$ is the weighting function for pressure. We represent these weighting functions discretely on a per-element-basis as

$$ w_{i} = \sum_{a=1}^{n_{en}} N_{a}^{w}w_{ai}, $$

$$ q = \sum_{a=1}^{n_{en}} N_{a}^{q}q_{a}, $$

where $N_{a}^{w}$ and $N_{a}^{q}$ are the nodal shape (basis) functions for the velocity and pressure spaces, respectively, and $w_{ai}$ and $q_{a}$ are the associated arbitrary nodal coefficients. Similarly, the trial functions are represented by

$$ u_{i} = \sum_{a=1}^{n_{en}} N_{a}^{w}u_{ai}, $$

$$ p = \sum_{a=1}^{n_{en}} N_{a}^{q}p_{a}. $$

We then multiply the Navier-Stokes-Brinkman equations by $\boldsymbol{w}$ and $q$, respectively, and integrate by parts to obtain the standard Galerkin momentum and continuity weak forms [7] [5],

$$ \int_{\Omega} \rho w_{i}\frac{du_{i}}{dt} d\Omega + \int_{\Omega} \rho w_{i}u_{k}u_{i, k} d\Omega + \int_{\Omega} w_{i, j}\sigma_{ij} d\Omega + \int_{\Omega} \frac{\mu}{K}w_{i}u_{i} d\Omega - \int_{\Omega} w_{i}\rho b_{i} d\Omega - \int_{\Gamma_{h}} w_{i}h_{i} d\Gamma = 0, $$

$$ \int_{\Omega} qu_{i,i} d\Omega = 0. $$

These two equations can be added together to obtain

$$ \int_{\Omega} qu_{i,i} d\Omega + \int_{\Omega} \rho w_{i}\frac{du_{i}}{dt} d\Omega + \int_{\Omega} \rho w_{i}u_{k}u_{i, k} d\Omega + \int_{\Omega} w_{i, j}\sigma_{ij} d\Omega + \int_{\Omega} \frac{\mu}{K}w_{i}u_{i} d\Omega - \int_{\Omega} w_{i}\rho b_{i} d\Omega - \int_{\Gamma_{h}} w_{i}h_{i} d\Gamma = 0. $$

Stabilized weak form

The standard weak form is generally not stable. Additional terms must be added to stabilize it. We will apply the residual-based variational multiscale (RBVMS / VMS) method for stabilization [1] [5].

In VMS, the velocity and pressure terms are separated into coarse-scale and fine-scale components, such that

$$ \boldsymbol{u} = \boldsymbol{u}^{h} + \boldsymbol{u}', \ $$

$$ p = p^{h} + p', $$

where the $h$-superscript designates the coarse-scale components and the $'$-superscript denotes the fine-scale components. The fine-scale terms are defined as

$$ \boldsymbol{u}' = -\frac{\tau_{SUPS}}{\rho}\boldsymbol{r}_{M}\left(\boldsymbol{u}^{h}, p^{h}\right), $$

$$ p' = -\rho\nu_{LSIC}r_{C}\left(\boldsymbol{u}^{h}\right), $$

where the PDE residuals are

$$ \boldsymbol{r}_{M}\left(\boldsymbol{u}^{h}, p^{h}\right) = \rho\left(\frac{d\boldsymbol{u}^{h}}{dt} + \boldsymbol{u}^{h} \cdot \boldsymbol{\nabla} \boldsymbol{u}^{h} - \boldsymbol{b}\right) - \boldsymbol{\nabla} \cdot \boldsymbol{\sigma}^{h} + \frac{\mu}{K}\boldsymbol{u}^{h}, $$

$$ r_{C}\left(\boldsymbol{u}^{h}\right) = \boldsymbol{\nabla} \cdot \boldsymbol{u}^{h}. $$

The stabilization parameters are defined as

$$ \tau_{SUPS} = \tau_{M} = \left(\frac{4}{\Delta t^{2}} + \boldsymbol{u}^{h} \cdot \boldsymbol{G}\boldsymbol{u}^{h} + C_{1}\nu^{2}\boldsymbol{G}:\boldsymbol{G} + \left(\frac{\nu}{K}\right)^{2}\right)^{-1/2}, $$

$$ \nu_{LSIC} = \tau_{C} = \left(\tau_{SUPS} \text{tr}\boldsymbol{G} \right)^{-1}, $$

where $\boldsymbol{G}$ is the element metric tensor and $\text{tr}\boldsymbol{G}$ is the trace of the metric tensor [1].

Using standard Galerkin momentum and continuity weak forms, and removing the $h$-superscript from the coarse-scale components for notational simplicity (i.e., $\boldsymbol{u}^{h} \rightarrow \boldsymbol{u}$ and $p^{h} \rightarrow p$), we obtain

$$ \int_{\Omega} qu_{i,i} d\Omega + \int_{\Omega} \rho w_{i}\frac{du_{i}}{dt} d\Omega + \int_{\Omega} \rho w_{i}u_{k}u_{i, k} d\Omega + \int_{\Omega} w_{i, j}\sigma_{ij} d\Omega + \int_{\Omega} \frac{\mu}{K}w_{i}u_{i} d\Omega - \int_{\Omega} w_{i}\rho b_{i} d\Omega - \int_{\Gamma_{h}} w_{i}h_{i} d\Gamma + \int_{\Omega} \tau_{SUPS}\left(\frac{q_{,i}}{\rho} + w_{i,k}u_{k}\right)r_{Mi} d\Omega + \int_{\Omega} \rho \nu_{LSIC}r_{C}w_{i,i} d\Omega - \int_{\Omega} w_{i}\tau_{SUPS}r_{Mk}u_{i,k} d\Omega - \int_{\Omega} w_{i,k}\frac{\tau_{SUPS}^{2}}{\rho}r_{Mi}r_{Mk} d\Omega - \int_{\Omega} \frac{\nu}{K}w_{i}\tau_{SUPS}r_{Mi} d\Omega = 0. $$

This is the VMS-stabilized weak form for the Navier-Stokes-Brinkman equations [7] [5] [6]. The first seven terms on the left-hand side correspond to the standard Galerkin weak form. The last five terms are the stabilization terms obtained via VMS. In deriving this equation, we used the continuity equation to obtain $w_{i}u_{k}u_{i,k} = w_{i}\left(u_{k}u_{i}\right)_{,k}$. We also applied the following assumptions [5],

We then add an additional stabilization term, $\int_{\Omega} \frac{\bar{\tau}\tau_{SUPS}^{2}}{\rho} w_{i,k}r_{Mk}r_{Mj}u_{i,j} d\Omega$ [3] [4] [2], to obtain

$$ \int_{\Omega} qu_{i,i} d\Omega + \int_{\Omega} \rho w_{i}\frac{du_{i}}{dt} d\Omega + \int_{\Omega} \rho w_{i}u_{k}u_{i, k} d\Omega + \int_{\Omega} w_{i, j}\sigma_{ij} d\Omega + \int_{\Omega} \frac{\mu}{K}w_{i}u_{i} d\Omega - \int_{\Omega} w_{i}\rho b_{i} d\Omega - \int_{\Gamma_{h}} w_{i}h_{i} d\Gamma + \int_{\Omega} \tau_{SUPS}\left(\frac{q_{,i}}{\rho} + w_{i,k}u_{k}\right)r_{Mi} d\Omega + \int_{\Omega} \rho \nu_{LSIC}r_{C}w_{i,i} d\Omega - \int_{\Omega} w_{i}\tau_{SUPS}r_{Mk}u_{i,k} d\Omega - \int_{\Omega} w_{i,k}\frac{\tau_{SUPS}^{2}}{\rho}r_{Mi}r_{Mk} d\Omega - \int_{\Omega} \frac{\nu}{K}w_{i}\tau_{SUPS}r_{Mi} d\Omega + \int_{\Omega} \frac{\bar{\tau}\tau_{SUPS}^{2}}{\rho} w_{i,k}r_{Mk}r_{Mj}u_{i,j} d\Omega = 0. $$

This equation is the full stabilized weak form used in fluid.cpp in svMultiPhysics.

Residuals

We will temporally discretize the stabilized weak form using the generalized - $\alpha$ method. The resulting nonlinear equation will be linearized and solved iteratively using the Newton-Raphson (Newton) method.

To compute the residuals for each element in the mesh, we separate the stabilized weak form into momentum and continuity components,

$$ \int_{\Omega} \rho w_{i}\frac{du_{i}}{dt} d\Omega + \int_{\Omega} \rho w_{i}u_{k}u_{i, k} d\Omega + \int_{\Omega} w_{i, j}\sigma_{ij} d\Omega + \int_{\Omega} \frac{\mu}{K}w_{i}u_{i} d\Omega - \int_{\Omega} w_{i}\rho b_{i} d\Omega - \int_{\Gamma_{h}} w_{i}h_{i} d\Gamma + \int_{\Omega} \tau_{SUPS}w_{i,k}u_{k}r_{Mi} d\Omega + \int_{\Omega} \rho \nu_{LSIC}r_{C}w_{i,i} d\Omega - \int_{\Omega} w_{i}\tau_{SUPS}r_{Mk}u_{i,k} d\Omega - \int_{\Omega} w_{i,k}\frac{\tau_{SUPS}^{2}}{\rho}r_{Mi}r_{Mk} d\Omega - \int_{\Omega} \frac{\nu}{K}w_{i}\tau_{SUPS}r_{Mi} d\Omega + \int_{\Omega} \frac{\bar{\tau}\tau_{SUPS}^{2}}{\rho} w_{i,k}r_{Mk}r_{Mj}u_{i,j} d\Omega = 0, $$

$$ \int_{\Omega} qu_{i,i} d\Omega + \int_{\Omega} \tau_{SUPS}\frac{q_{,i}}{\rho}r_{Mi} d\Omega = 0. $$

Then, by plugging weighting functions into these equations and holding the results true for any arbitrary $w_{ai}$ and $q_{a}$, we obtain the momentum and continuity residuals,

$$ R_{ai}^{m} = \int_{\Omega} \rho N_{a}^{w}\frac{du_{i}}{dt} d\Omega + \int_{\Omega} \rho N_{a}^{w}u_{k}u_{i, k} d\Omega - \int_{\Omega} pN_{a, i}^{w} d\Omega + \int_{\Omega} N_{a, j}^{w}2\mu\epsilon_{ij} d\Omega + \int_{\Omega} \frac{\mu}{K}N_{a}^{w}u_{i} d\Omega - \int_{\Omega} N_{a}^{w} \rho b_{i} d\Omega + \int_{\Omega} \tau_{SUPS}N_{a, k}^{w}u_{k}r_{Mi} d\Omega + \int_{\Omega} \rho \nu_{LSIC}r_{C}N_{a, i}^{w} d\Omega - \int_{\Omega} N_{a}^{w}\tau_{SUPS}r_{Mk}u_{i,k} d\Omega - \int_{\Omega} N_{a, k}^{w}\frac{\tau_{SUPS}^{2}}{\rho}r_{Mi}r_{Mk} d\Omega - \int_{\Omega} \frac{\nu}{K}\tau_{SUPS}N_{a}^{w}r_{Mi} d\Omega + \int_{\Omega} \frac{\bar{\tau}\tau_{SUPS}^{2}}{\rho} N_{a, k}^{w}r_{Mk}r_{Mj}u_{i,j} d\Omega, $$

$$ R_{a}^{c} = \int_{\Omega} N_{a}^{q}u_{i,i} d\Omega + \int_{\Omega} \tau_{SUPS}\frac{N_{a, i}^{q}}{\rho}r_{Mi} d\Omega, $$

where, for the $a^{\text{th}}$ node in a given element, $R_{ai}^{m}$ is the momentum residual in the $i^{\text{th}}$ direction and $R_{a}^{c}$ is continuity residual. The full residual vector, as used in the generalized - $\alpha$ method, is $\boldsymbol{R} = \left[R_{ai}^{m}, R_{a}^{c}\right]^{T}$. $R_{ai}^{m}$ and $R_{a}^{c}$ are coded in the fluid_2d_m/fluid_3d_m and fluid_2d_c/fluid_3d_c functions, respectively, in fluid.cpp in svMultiPhysics.

Tangent matrices

To compute the elemental tangent matrices, as used in the Newton iterations in the generalized - $\alpha$ method, we plug trial functions into the residuals. We then differentiate the resulting equations with respect to $\frac{du_{n+1}}{dt}$ and $\frac{dp_{n+1}}{dt}$. This yields the tangent matrix,

$$ \boldsymbol{J} = \begin{bmatrix} \boldsymbol{K} & \boldsymbol{G} \ \cr \boldsymbol{D} & \boldsymbol{L} \end{bmatrix} = \begin{bmatrix} \left[K_{ab}^{ij}\right] & \left[G_{ac}^{i}\right] \ \cr \left[D_{ab}^{j}\right] & \left[L_{ac}\right] \end{bmatrix} , $$

where

$$ K_{ab}^{ij} = \frac{\partial R_{ai}^{m}}{\partial \dot{u}_{n+1,bj}}, $$

$$ G_{ac}^{i} = \frac{\partial R_{ai}^{m}}{\partial \dot{p}_{n+1,c}}, $$

$$ D_{ab}^{j} = \frac{\partial R_{a}^{c}}{\partial \dot{u}_{n+1,bj}}, $$

$$ L_{ac} = \frac{\partial R_{a}^{c}}{\partial \dot{p}_{n+1,c}}. $$

In fluid.cpp of svMultiPhysics, the following inconsistent tangent matrices are used,

$$ K_{ab}^{ij} = \alpha_{m} A_{ab}^{ij} + \alpha_{f}\gamma\Delta t B_{ab}^{ij} $$

$$ G_{ac}^{i} = \alpha_{f}\gamma\Delta t \left(-\int_{\Omega} N_{c}^{q}N_{a, i}^{w} d\Omega + \int_{\Omega} \tau_{SUPS} N_{a, g}^{w} u_{g} N_{c, i}^{q} d\Omega - \int_{\Omega} N_{a, k}^{w} \frac{\tau_{SUPS}^{2}}{\rho} N_{c, i}^{q} r_{Mk} d\Omega \right), $$

$$ D_{ab}^{j} = \alpha_{f}\gamma\Delta t \left(\int_{\Omega} N_{a}^{q}N_{b, j}^{w} d\Omega - \int_{\Omega} \tau_{SUPS}\frac{N_{a, i}^{q}}{\rho}\left(-\frac{\alpha_{m}}{\alpha_{f}\gamma\Delta t}\rho N_{b}^{w}\delta_{ij} - \frac{\partial r_{Mi}}{\partial u_{n+\alpha_f,bj}}\right) d\Omega\right), $$

$$ L_{ac} = \alpha_{f}\gamma\Delta t \int_{\Omega} \tau_{SUPS}\frac{N_{a, i}^{q}}{\rho}N_{c, i}^{q} d\Omega, $$

where

$$ A_{ab}^{ij} = \int_{\Omega} \left( \rho N_{a}^{w}N_{b}^{w} \delta_{ij} + \tau_{SUPS} N_{a,g}^{w} u_{g} \rho N_{b}^{w} \delta_{ij} - N_{a,k}^{w} \tau_{SUPS}^{2} N_{b}^{w} \delta_{ij} r_{Mk} \right) d\Omega , $$

$$ B_{ab}^{ij} = \int_{\Omega} \left( \rho N_{a}^{w} u_{k} N_{b, k}^{w} \delta_{ij} + N_{a, l}^{w} \mu N_{b, l}^{w} \delta_{ij} + N_{a, j}^{w} \mu N_{b, i}^{w} + \frac{\mu}{K} N_{a}^{w} N_{b}^{w} \delta_{ij} + \tau_{SUPS} N_{a,g}^{w} u_{g} \frac{\partial r_{Mi}}{\partial u_{n+\alpha_f,bj}} + \rho \nu_{LSIC} N_{b,j}^{w} N_{a,i}^{w} - N_{a}^{w} \tau_{SUPS} N_{b,k}^{w} \delta_{ij} r_{Mk} - N_{a,k}^{w} \frac{\tau_{SUPS}^{2}}{\rho} \frac{\partial r_{Mi}}{\partial u_{n+\alpha_f,bj}} r_{Mk} + \frac{4}{\gamma} \frac{\partial \mu}{\partial \gamma} \epsilon_{jk} N_{b,k}^{w} \epsilon_{il} N_{a,l}^{w} + \frac{\bar{\tau}\tau_{SUPS}^{2}}{\rho} N_{a,k}^{w} N_{b,z}^{w} r_{Mk} r_{Mz} \delta_{ij} \right) d\Omega , $$

and

$$ \frac{\partial r_{Mi}}{\partial u_{n+\alpha_f,bj}} = \left(\rho u_{k} N_{b,k}^{w} - \mu N_{b,kk}^{w} + \frac{\mu}{K} N_{b}^{w} - \frac{\partial \mu}{\partial x_{k}} N_{b,k}^{w} \right)\delta_{ij} - \frac{2}{\gamma} \frac{\partial \mu}{\partial \gamma} \epsilon_{il} N_{b,l}^{w} u_{j, kk} - \frac{\partial \mu}{\partial x_{j}} N_{b,i}^{w}. $$

These inconsistent tangent matrices were derived by using these assumptions:

$K_{ab}^{ij}$ and $G_{ac}^{i}$ are coded in the fluid_2d_m/fluid_3d_m functions, while $D_{ab}^{j}$ and $L_{ac}$ are coded in the fluid_2d_c/fluid_3d_c functions.

References

[1] Y. Bazilevs, V.M. Calo, J.A. Cottrell, T.J.R. Hughes, A. Reali, and G. Scovazzi (2007). Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows. Computer Methods in Applied Mechanics and Engineering, 197, 173-201, https://doi.org/10.1016/j.cma.2007.07.016

[2] C.A. Taylor, T.J.R. Hughes, and C.K. Zarins (1998). Finite element modeling of blood flow in arteries. Computer Methods in Applied Mechanics and Engineering, 158, 155-196, https://doi.org/10.1016/S0045-7825(98)80008-X

[3] M. Esmaily-Moghadam (2014). Development of multiscale modeling methods for clinical decision making in single ventricle heart patients. University of California, San Diego

[4] C.H. Whiting and K.E. Jansen (2001). A stabilized finite element method for the incompressible Navier–Stokes equations using a hierarchical basis. International Journal for Numerical Methods in Fluids, 35, 93-116, https://doi.org/10.1002/1097-0363(20010115)35:1<93::AID-FLD85>3.0.CO;2-G

[5] Y. Bazilevs, K. Takizawa, T.E. Tezduyar (2013). Computational Fluid–Structure Interaction: Methods and Applications. John Wiley & Sons, Ltd.

[6] F.M. Gerosa and A.L. Marsden (2024). A mechanically consistent unified formulation for fluid-porous-structure-contact interaction. Computer Methods in Applied Mechanics and Engineering, 425, 116942, https://doi.org/10.1016/j.cma.2024.116942

[7] J. Fuchsberger, P. Aigner, S. Niederer, G. Plank, H. Schima, G. Haase, and E. Karabelas (2022). On the incorporation of obstacles in a fluid flow problem using a Navier–Stokes–Brinkman penalization approach. Journal of Computational Science, 57, 101506, https://doi.org/10.1016/j.jocs.2021.101506

[8] T.J.R. Hughes (2003). CThe Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover Publications

Solver Parameter Input File

The svMultiPhysics solver reads simulation parameters from a text file in Extensible Markup Language (XML) format. The XML file structurally organizes the svMultiPhysics solver parameters and input data (e.g. mesh files) needed to set up and execute a simulation.

What is XML

The XML format is made up of tags, elements and attributes. A tag begins with < and ends with >. There are two types of tag

• start-tag, such as <section>
• end-tag , such as </section>

An XML element is everything from (including) the element's start tag to (including) the element's end tag.

<Tolerance> 0.001 </Tolerance>

There are three types of elements

• text
• attributes - name-value pair within a start-tag: name="value" "
• other elements - elements nested within other elements

The following XML is an example of these different element types

<Add_mesh name="fluid">


 <Mesh_file_path> fluid_mesh.vtu 
 </Mesh_file_path>




 <Add_face name="inlet">

  <Face_file_path> inlet.vtp
</Face_file_path>


 </Add_face>





 <Domain> 1 
</Domain>







 <Mesh_scale_factor> 0.1 
</Mesh_scale_facto>



</Add_mesh>

<?xml version="1.0" encoding="UTF-8" ?>


<svMultiPhysicsFile version="0.1">

</svMultiPhysicsFile>

<GeneralSimulationParameters>
[general parameters]
</GeneralSimulationParameters>


<Add_mesh name=mesh_name_1>
  <Add_face name=mesh1_face_name_1>
  <Add_face name=mesh1_face_name_2>
  ...
  <Add_face name=mesh1_face_name_n>
</Add_mesh>

<Add_mesh name=mesh_name_2>
  <Add_face name=mesh2_face_name_1>
  <Add_face name=mesh2_face_name_2>
  ...
  <Add_face name=mesh2_face_name_n>
</Add_mesh>
...
<Add_mesh name=mesh_name_N>
  <Add_face name=meshN_face_name_1>
  <Add_face name=meshN_face_name_2>
  ...
  <Add_face name=meshN_face_name_M>
</Add_mesh>


<Add_equation type=eq_1>
  [equation parameters]

  <Domain id=0>
  <Domain id=1>
  ...
  <Domain id=N>

  <LS type=ls_type>
  [linear solver parameters]
  </LS>

  <Output type=output_type>
  [solver output parameters]
  </Output>

</Add_equation>

<Spectral_radius_of_infinite_time_step> real  [0.5] 
</Spectral_radius_of_infinite_time_step>

<Save_results_in_folder> string [none] 
</Save_results_in_folder>

<Number_of_time_steps> integer 
</Number_of_time_steps>

General Simulation Parameters

• spatial dimension
• time step control
• results output

Parameters

Mesh Section

A Mesh Section is organized as follows

The <Add_mesh name=mesh_name> parameter adds a volumetric mesh named mesh_name to the simulation. Multiple <Add_mesh name=mesh_name> parameters can be given within a solver parameters input file.

The <Add_face name=mesh_face_name> parameters subsection under the <Add_mesh> parameter associates the surface mesh named mesh_face_name with the volumetric mesh. A surface mesh, referred to as a face, is used to apply boundary conditions.

Mesh Parameters

Add_Face Parameters

Contact Section

A Contact Section is organized as follows

The Contact keyword is used to define the parameters used in a contact computation. The model attribute contact_model is the name of a method used for the contact computation. The "penalty" model is currently the only model supported.

Contact Parameters

Precomputed Solution Section

A Precomputed Solution Section is organized as follows

The Precomputed_solution keyword is used to define the parameters to for reading in a precomputed solution data and using it in a simulation.

Precomputed Solution Parameters

Projection Section

A Projection Section is organized as follows

The Projection keyword is used to create a projection. The name attribute project_from_face_name is the name of a face defined in the Mesh Section.

Projection Parameters

Equation Section

• physics - fluid, structure, electrophysiology, etc.
• domains

The Equation Section is organized as follows

The Add_equation keyword adds an equation of type equation_type to the simulation. Multiple Add_equation keywords can be given within a solver parameters input file.

The value of equation_type is a string identifying the equation type

• "advection_diffusion" - unsteady advection-diffusion
• "cardiac_electro_physiology" - solves the mono-domain model of cardiac electrophysiology
• "coupled_momentum" - coupled momentum method for fluid-structure interaction
• "fluid" - Navier-Stokes equations for unsteady viscous incompressible fluid flow
• "fluid-solid-interaction" - fluid-solid interaction using the arbitrary Lagrangian-Eulerian (ALE) formulation
• "linear_elasticity" - linear elastodynamics equation
• "mesh" - solves a modified linear elasticity equation for mesh motion in an fluid-solid interaction simulation
• "shell" - nonlinear thin shell mechanics (Kirchhoff-Love theory)
• "solid_heat" - unsteady diffusion equation
• "stokes" - unsteady Stokes equations
• "structural" - nonlinear elastodynamics equation
• "structural_velocity_pressure" - nonlinear elastodynamics equation using mixed VMS-stabilized formulation

Equation Parameters

ECG Leads Subsection

The ECG Leads Subsection is organized as follows

The ECGLeads keyword defines a subsection for ECG Leads parameters.

ECG Leads Parameters

Variable wall properties Subsection

The Variable wall properties Subsection is organized as follows

The Variable_wall_properties keyword defines a subsection for variable wall properties parameters. The mesh_name attribute is the mesh name to assign the variable wall properties.

Variable wall properties Parameters

Viscosity Subsection

The Viscosity Subsection is organized as follows

The Viscosity keyword defines a subsection for viscosity parameters.

For fluid, CMM, and stokes equations, the value of viscosity_model can be

• "newtonian" - Newtonian viscosity model
• "carreau-yasuda" - Carreau-Yasuda viscosity model
• "cassons" - Cassons viscosity model

Newtonian

Carreau-Yasuda

Cassons

For struct and ustruct equations, the value of viscosity_model can be

• "Newtonian" - Newtonian viscosity model
• "Potential" - Pseudo-potential viscosity model

Newtonian

Potential

Domain Subsection

The Domain Subsection is organized as follows

The Domain keyword adds a domain to the enclosing equation. domain_id is an integer ID defined in Mesh Section for a finite element mesh.

Domain Parameters by Equation

Advection Diffusion Domain Parameters

 Conductivity 
 Source_term 

 Anisotropic_conductivity 
 Electrophysiology_model 
 Myocardial_zone 
 ODE_solver 
 Stimulus 
 Time_step_for_integration 

 Elasticity_modulus 
 Fluid_density 
 Solid_density
 Backflow_stabilization_coefficient 
 Poisson_ratio 
 Shell_thickness 

 Backflow_stabilization_coefficient 
 Density 
 Force_x 
 Force_y 
 Force_z 
 Inverse_darcy_permeability 

 Density 
 Elasticity_modulus 
 Force_x 
 Force_y 
 Force_z 
 Poisson_ratio 

 Poisson_ratio 

 Density 
 Elasticity_modulus 
 Force_x 
 Force_y 
 Force_z 
 Poisson_ratio 

 Conductivity 
 Density 
 Source_term 

 Force_x 
 Force_y 
 Force_z 
 Momentum_stabilization_coefficient 

 Density 
 Elasticity_modulus 
 Force_x 
 Force_y 
 Force_z 
 Poisson_ratio 

 Density 
 Elasticity_modulus 
 Force_x 
 Force_y 
 Force_z 
 Momentum_stabilization_coefficient 
 Poisson_ratio 

The Stimulus keyword adds stimulus settings to the enclosing domain. The stimulus_type attribute can be

• "Istim" - the applied stimulus is a source current and not a voltage source

Parameters

Linear Solver Subsection

The Linear Solver Subsection is organized as follows

The LS keyword adds a linear solver to the enclosing equation.

The value of linear_solver_type can be

Linear Solver Parameters

Linear Algebra Parameters

The value of linear_algebra_type can be

• "fsils" - Custom numerical linear algebra routines included in the svMultiPhysics implementation
• "petsc" - Portable, Extensible Toolkit for Scientific Computation (PETSc)
• "trilinos" - The Trilinos numerical linear algebra package

svZeroDSolver Interface Subsection

The SimVascular svZeroDSolver simulates bulk cardiovascular flow rates and pressures using an arbitrary zero-dimensional (0D) lumped parameter model (LPM) of a discrete network of components analogous to electrical circuits. It provides an Application Programming Interface (API) that allows it to communicate and interact with external software applications directly using function calls to programmatically define custom inflow and outflow boundary conditions for a CFD simulation. The svMultiPhysics solver can directly access the svZeroDSolver API by loading the svZeroDSolver as a shared (dynamic) library available after installing the svZeroDSolver.

The svZeroDSolver Interface Subsection is organized as follows

The svZeroDSolver_interface keyword activates coupling for boundary conditions with the <Time_dependence> Coupled </Time_dependence> keyword for the enclosing equation coupled boundary condition.

The value of svzerodsolver_interface_name can be any string name and is not currently used.

svZeroDSolver Interface Parameters

Output Subsection

The Output Subsection is organized as a list of quantity names with a boolean

The Output keyword enables output of quantities for the enclosing equation. The value of output_type can be

• "Boundary_integral" - Quantities computed as a flux through boundary surfaces (faces) (e.g. velocity flux, energy flux)
• "Spatial" - Quantities written to VTK VTU files
• "Volume_integral" - Volume-averaged quantities

Quantity names are boolean parameters used to enable quantities for output by setting their value to true. The quantities available for output depends on the equation being solved. Below is a list of all quantity names available for output for the corresponding equation

• advection_diffusion - Temperature, Heat_flux
• cardiac_electro_physiology - Action_potential
• coupled_momentum - Velocity, Pressure, Displacement,WSS, Vorticity, Vortex, Energy_flux, Strain_invariants, Acceleration, Traction, Viscosity, Divergence
• fluid - Velocity, Pressure, WSS, Traction, Vorticity, Vortex, Energy_flux, Strain_invariants, Acceleration, Viscosity, Divergence
• fluid-solid-interaction - Velocity, Pressure, Displacement, VonMises_stress, WSS, Traction, Vorticity, Vortex, Energy_flux, Strain_invariants, Viscosity, Absolute_velocity, Stress, Strain, Cauchy_stress, Jacobian, Deg_grad, Area/Volume, Fiber_direction, Fiber_alignment, Divergence, Acceleration
• linear_elasticity - Displacement, VonMises_stress, Stress, Strain, Jacobian, Area/Volume, Velocity, Acceleration
• mesh - Displacement, Velocity, Acceleration
• shell - Displacement, Stress, Strain (Green-Lagrange in-plane), Def_grad, Area, Jacobian, Velocity, 1st invariant of Cauchy-Green strain (in-plane), Cauchy-Green strain tensor
• solid_heat - Temperature, Heat_flux
• stokes - Velocity, Pressure, WSS, Vorticity, Traction, Strain_invariants, Viscosity, Divergence
• structural - Displacement, VonMises_stress, Stress, Cauchy_stress, Strain, Jacobian, Def_grad, Area/Volume, Fiber_direction, Fiber_alignment, Velocity, Acceleration
• structural_velocity_pressure - Displacement, VonMises_stress, Stress, Cauchy_stress, Strain, Jacobian, Def_grad, Area/Volume, Fiber_direction, Fiber_alignment, Velocity, Pressure, Acceleration, Divergence

Example

Boundary Condition Parameters

The Boundary Condition Parameters section is organized as follows

The <Add_BC name=face_name> parameter adds a boundary condition to the enclosing equation.

Boundary condition types

Temporal and spatial values for boundary conditions

General Parameters

Dirichlet Boundary Condition Parameters

Neumann Boundary Condition Parameters

Robin Boundary Condition Parameters

Boundary Condition Data

• Fourier coefficients
• Spatial distribution
• Temporal distribution
• Temporal and spatial distribution
• User-defined flow profile

• ASCII text file usually with a .dat file extension.
• VTK format file

Fourier coefficients Data

This file is set using the solver input file Fourier_coefficients_file_path keyword.

The file format is

InitialTime Period
InitialValue_1 InitialValue_2 ... InitialValue_NDIM 
TimeDerivativeLinearPart1 TimeDerivativeLinearPart2 ... TimeDerivativeLinearPart_NDIM
NumberOfFourierModes
RealPart_1 ImaginaryPart_1
...

where NDIM is the number of space dimensions used in the simulation.

Spatial Distribution Data

• Neumann boundary condition
• Traction boundary condition
• Coupled Momentum boundary condition

Neumann and Traction boundary conditions

• Pressure - Neumann boundary condition
• Traction - Traction boundary condition

A spatial data file is set using the solver input file Spatial_values_file keyword.

Coupled momentum boundary condition

• Stress - Initialize stresses
• Displacement - Initialize displacements

These spatial data files are set using the solver input file Initial_displacements_file_path and Prestress_file_path and keywords.

Temporal Distribution Data

Temporal data will be interpolated internally using Fourier series interpolation to obtain periodic data at intermediate time values.

This file is set using the solver input file Temporal_values_file_path keyword.

The file format is

NumberOfTimePoints  NumberOfFourierModes  
Time_1  Value_1  
Time_2  Value_2  
...
Time_NumberOfTimePoints  Value_NumberOfTimePoints  

Example: Inlet fluid flow rate

This file is set using the solver input file Temporal_and_spatial_values_file_path keyword.

The file format is

NumberDegreesOfFreedom NumberOfTimePoints  NumberOfNodes
Time_1  
Time_2 
...
Time_NumberOfTimePoints  
NodeID_1
Time_1  Value_1,1 Value_1,2  ... Value_1,NumberDegreesOfFreedom
Time_2  Value_2,1 Value_2,2  ... Value_2,NumberDegreesOfFreedom
....
NodeID_2
Time_1  Valu_1 Value_2  ... Value_NumberDegreesOfFreedom
...
NodeID_NumberOfNodes
...

Example: Volumetric body force

This file is set using the solver input file Spatial_profile_file_path keyword.

The file format is

NodeID_1 Value_1
NodeID_2 Value_2
...
NodeID_NumberOfFaceNodes Value_NumberOfFaceNodes 

where NumberOfFaceNodes is the number of nodes in the boundary surface.

Fiber and Sheet Direction Data

The data file is set using the solver input file Fiber_direction_file_path keyword.

Initial Condition Data

Displacement, Pressure, Stress, and Velocity Data

The data files are set using the solver input file Initial_displacements_file_path, Initial_pressures_file_path and Prestress_file_path and Initial_velocities_file_path keywords.

State Variables

The data file is set using the solver input file Simulation_initialization_file_path keyword.

Finite Element Mesh Data

Finite Element Volume Mesh

• number of nodes - NumberOfPoints
• nodal coordinates - DataArray type="Float32" Name="Points"
• number of elements - NumberOfCells
• element connectivity - DataArray type="Int64" Name="connectivity"
• user-defined arrays - CellData (element) or PointData (node) data specific to a SimVascular application
• GlobalElementID - CellData
• GlobalNodeID - PointData
• ModelRegionID - CellData

A finite element volume mesh file is set using the solver input file Mesh_file_path keyword.

Example: A finite element volume mesh comprising 8561 nodes and 45723 elements

<?xml version="1.0"?>
<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
  <UnstructuredGrid>
    <Piece NumberOfPoints="8561"                 NumberOfCells="45723"               >
      <PointData>
        <DataArray type="Int32" Name="GlobalNodeID" format="appended" RangeMin="1"                    RangeMax="8561"                 offset="0"                   />
      </PointData>
      <CellData Scalars="ModelRegionID">
        <DataArray type="Int32" Name="ModelRegionID" format="appended" RangeMin="1"                    RangeMax="1"                    offset="15964"               />
        <DataArray type="Int32" Name="GlobalElementID" format="appended" RangeMin="1"                    RangeMax="45723"                offset="16444"               />
      </CellData>
      <Points>
        <DataArray type="Float32" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0.12837002719"        RangeMax="30.066622254"         offset="100980"              />
      </Points>
      <Cells>
        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="227180"              />
        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="788628"              />
        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="868152"              />
      </Cells>
    </Piece>
  </UnstructuredGrid>
[COMPRESSED DATA]

The nodes and elements in the finite element mesh are identified by a unique ID in the GlobalNodeID and GlobalElementID arrays, respectively.

Finite Element Boundary Surface Mesh

• number of nodes - NumberOfPoints
• nodal coordinates - DataArray type="Float32" Name="Points"
• number of face elements - NumberOfCells
• element face connectivity - DataArray type="Int64" Name="connectivity"
• user-defined arrays - CellData (element) or PointData (node) data specific to a SimVascular application
• GlobalElementID - CellData
• GlobalNodeID - PointData
• ModelFaceID - CellData

A finite element boundary surface mesh file is set using the solver input file Face_file_path keyword.

Example: A finite element boundary surface mesh comprising 2647 nodes and 5290 elements

<?xml version="1.0"?>
<VTKFile type="PolyData" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
  <PolyData>
    <Piece NumberOfPoints="2647"                 NumberOfVerts="0"                    NumberOfLines="0"                    NumberOfStrips="0"                    NumberOfPolys="5290"                >
      <PointData Scalars="GlobalNodeID">
        <DataArray type="Int32" Name="GlobalNodeID" format="appended" RangeMin="1"                    RangeMax="2647"                 offset="0"                   />
      </PointData>
      <CellData Scalars="ModelFaceID">
        <DataArray type="Int32" Name="GlobalElementID" format="appended" RangeMin="2"                    RangeMax="45702"                offset="4964"                />
        <DataArray type="Int32" Name="ModelFaceID" format="appended" RangeMin="1"                    RangeMax="3"                    offset="23712"               />
      </CellData>
      <Points>
        <DataArray type="Float32" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0.12837002719"        RangeMax="30.066622254"         offset="23932"               >
        </DataArray>
      </Points>
      <Verts>
        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="62272"               />
        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="62288"               />
      </Verts>
      <Lines>
        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="62304"               />
        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="62320"               />
      </Lines>
      <Strips>
        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="62336"               />
        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="62352"               />
      </Strips>
      <Polys>
        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="62368"               />
        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="99592"               />
      </Polys>
    </Piece>
  </PolyData>
  <AppendedData encoding="base64">
[COMPRESSED DATA]

The node IDs contained in the GlobalNodeID array identify the nodes stored in the boundary surface VTP with the nodes in the finite element volume mesh.

The element IDs contained in the GlobalElementID array identify the element in the finite element volume mesh that each polygon of the face belongs to.

Results Data

• Binary results restart format
• VTK results format

See the Simulation Results Output section for a description of the solver input file keywords used to control how results files are written.

Binary Results Restart File Format

Each restart file contains a header of seven integer (4-byte) values containing

1. Number of processors used to run the simulation
2. Number of equations
3. Number of meshes
4. Number of nodes
5. Number of unknowns in the coupled 0D-3D problem
6. Number of degrees of freedom
7. An error flag

These values are followed by an integer (4-byte) and two double (8-byte) values containing

1. Simulation time step
2. Simulation time
3. Wall clock time

Header values are checked when a simulation is restarted to ensure that the restart is valid for the the current simulation setup (e.g. same number of processors, same number of equations, etc.).

State variable data (e.g. velocity) is then written after the header. The results from all processors in a parallel simulation are written to the restart file..

VTK Results File Format

A separate VTK file is used to store results for a single time step. Results can be visualized using the ParaView visualization program.

Data from all processors in a parallel simulation are written to the VTK files.

SimVascular MultiPhysics Tool

The SimVascular SV MultiPhysics Tool is used to interactively create the mesh and solver XML input files needed to run an svMultiPhysics simulation.

The MultiPhysics Tool stores data in the Project/svMultiPhysics/Data Node directory. The MultiPhysics Tool stores values input via the GUI controls into the ProjectMultiPhysics/Data Node.multiphysicsjob file.

The following sections describe how to use the SimVascular MultiPhysics Tool GUI to set the mesh, boundary conditions, material properties and solver parameters needed for a simulation.

Creating Mesh Files

Creating an MultiPhysics Tool Instance

Type fluid_simulation into the job Name DialogBox. This names is used to identify the MultiPhysics Tool instance and to name the files and directories stored under the SimVascular project's MultiPhysics directory. Selecting OK creates an MultiPhysics Tool instance node under the SV Data Manager. Selecting this instance displays the MultiPhysics Tool panel.

The panel contains four buttons at the top used to primarily to create new svMultiPhysics simulations (jobs)

New job Button - Selecting to create a new MultiPhysics instance

Save job Button - Select to write the values of GUI controls into the SVPROJECT/MultiPhysicsJob Name.multiphysicsjob file.

Load job Button - Select to display a file browser used to select a .multiphysicsjob file. The values of GUI controls are then set from the values in this file.

Load mesh - Reads in the finite element mesh stored in the Project/Meshes directory. At startup SimVascular does not read in the finite element mesh so this must be initiated manually.

The panel contains four sub-panels used to input a specific category of data needed to create the files needed to run a svMultiPhysics solver simulation

• Domains
• Physics
• Simulation Parameters
• Run Simulation

A selecting a sub-panel name brings up the sub-panel's controls. The following sections describe how each of the sub-panels are used.

Physics Panel

Panel Layout

Usage

Add or remove equations - This GUI control is used to add equations to the simulation. Selecting the equation name in the left panel and selecting the > Button moves it to the right panel and adds that equation to the simulation. Use the < Button to remove that equation from the simulation.

When an equation is added the Buttons at the bottom of this control (Properties, Output, etc.) are used to set the parameters for the selected equation. The GUI controls for setting parameters for each of these are displayed directly beneath them.

Properties Button - Select to display the GUI controls used to set the physicals properties (e.g. density) needed by the equation.

Output Button - Select to display the GUI controls to add/remove names of the quantities for output.

BCs Button - Select to display the GUI controls to add/remove names of the quantities for output.

Advanced Button - Select to display a GUI panel to set advanced solver options.

Linear Solver Button - Select to display a GUI panel to set linear solver options.

Remesher Button - Select to display a GUI panel to set remeshing options.

GUI controls for setting parameters

Properties Button - Setting physical properties

The physical properties displayed in the panel depend of the selected equation.

Fluid equation properties

Density TextBox - Fluid density
Viscosity TextBox - Fluid viscosity

Linear elasticity equation properties

Density TextBox - Solid density
Elastic modulus TextBox - Elastic modulus
Poisson ratio TextBox - Poisson ratio. Value <= 0.5

Structure equation properties

Density TextBox - Solid density
Elastic modulus TextBox - Elastic modulus
Poisson ratio TextBox - Poisson ratio. Value <= 0.5
Constitutive model ComboBox - Constitutive model selected from

• stVK - St. Venant-Kirchhoff material model
• nHK - NeoHookean material model

FSI equation properties

Density TextBox - Fluid density
Viscosity TextBox - Fluid viscosity
Solid Density TextBox - Solid density
Poisson ratio TextBox - Poisson ratio. Value <= 0.5
Constitutive model ComboBox - Constitutive model selected from

• stVK - St. Venant-Kirchhoff material model
• nHK - NeoHookean material model

Output Button - Adding output quantities

The quantities output by default are displayed in the right sub-panel.

Selecting the output quantity name in the left panel and selecting the > Button moves it to the right panel and adds that quantity for output.

BCs Button - Adding Boundary Conditions

Add Button - Select to add a boundary condition. A Boundary Condition Setup popup window is displayed with a list of mesh boundary faces. Selecting a face name allows for setting the various parameters for a given type of boundary condition.

Boundary Condition Setup Parameters
Face list: ListTable - Lists the boundary surface face names defined for a domain. Selecting a face name adds a boundary condition for that face.

BC Type: ComboBox - Select boundary condition type from the list [ Neumann Dirichlet]

Directional BC CheckBox - If selected

Time dependence

Steady CheckBox - If selected

Value TextBox -

Unsteady CheckBox - If selected

Temporal values file FileBrowser - Select to identify the time values

Resistance - CheckBox - If selected

Value TextBox -

Coupled - CheckBox - If selected

General - CheckBox - If selected

Temporal and spatial values file FileBrowser - Select to identify the time values

Projection - CheckBox - If selected

Profiledependance

Zero out perimeter - CheckBox - If selected

Impose flux - CheckBox - If selected

Flat - CheckBox - If selected

Parabolic - CheckBox - If selected

User defined - CheckBox - If selected

Spatial profile file path FileBrowser - Select to identify the time values

Impose on state variable integral - CheckBox - If selected

Effective direction TextBox -

Reset Button -

Cancel Button -

OK Button -

Advanced Button - Setting advanced solver parameters

Coupled - CheckBox - If selected

Tolerance TextBox -

Residual dB reduction SpinBox -

Min iterations SpinBox -

Max iterations SpinBox -

Backflow stabilization SpinBox -

Reset Button -

Linear solverButton - Setting linear solver parameters

Coupled - CheckBox - If selected

Type ComboBox - Select the linear solver to use.

Max iterations TextBox -

Tolerance TextBox -

Krylov dimension TextBox -

Absolute Tolerance TextBox -

Preconditioner ComboBox - Select the linear solver preconditioner to use.

Simulation Parameters Panel

Panel Layout

Usage

Time control
 
Start from previous simulation CheckBox - If selected then start the simulation from 
the saved state of a previous simulation.

Number of time steps TextBox - The number of time steps to run the simulation.

Time step size TextBox - The simulation time step.

Restart files
 
Name TextBox - The name prefix of the files used to store simulation state data.

Save increment TextBox - Indicates how often to save restart files.

Save results after time step TextBox - Start saving simulation results after this time.

Save as VTK CheckBox - If selected then save simulation output to VTK files.

Name TextBox - The name prefix of the VTK files used to store output data.

Save increment TextBox - Indicates how often to save the VTK files.

Time-averaged results CheckBox - If selected then compute time-averaged results for the entire simulation using the using VTU files written during the simulation.

Advanced options CheckBox - If selected then enable setting advanced parameters.

Spectral radius of infinite time step TextBox - The spectral radius is used to compute parameters for the generalized alpha method. A value of 0.0 leads to an over-damped system while 1.0 leads to an undamped system.

Remeshing CheckBox - If selected then the finite element mesh for FSI simulations will be
automatically remeshed during a simulation.

Warning CheckBox - If selected then print out warning messages indicating unexpected
inputs or simulation behavior.

Verbose CheckBox - If selected print detailed messages when processing input.

Debug CheckBox - If selected the print additional information that may be used for debugging purposes.

Run Simulation Panel

Panel Layout

Usage

Create Input File Button - Select to create solver input files.

Number of processes Slider - Select the number of processors (cores) for MPI to use for a parallel simulation.

Run Button - Select to start a simulation.

Stop Button - Select to stop a simulation.

Fluid-Structure Interaction Modeling using svMultiPhysics

Cardiovascular applications that involve significant deformation are difficult to simulate and model. Some examples include vascular simulations where the wall dilates to more than 10% of the vessel diameter, cardiac simulations where the heart wall displaces as it contracts and expands, or coronary simulations where the coronary vessels displace with the motion of the heart.

One way to computationally model situations like these is to use a fluid-structure interaction (FSI) solver. In FSI, separate domains are defined for the fluid part and solid parts of the computational geometry. The respective equations governing fluid flow (typically Navier-Stokes) are solved in the fluid domain, while the equations governing solid mechanics are solved in the solid domain. The two domains then interact through their interface where the solution variables (displacements, velocities, pressures, stresses) are required to match. This interface acts as a coupled boundary condition for both domains, as the solution of one domain will affect the solution in the other, and vice-versa.

svMultiPhysics utilizes Arbitrary Lagrangian-Eulerian (ALE) method to perform FSI. As the name implies, this method is a combination of the Eulerian and Lagrangian descriptions of motion that is particularly well suited for FSI problems. In the Eulerian description, the computational mesh stays fixed and the motion is characterized as velocities flowing past the grid nodes. The Eulerian description is typical for most fluid dynamics problems without FSI. While in the Lagrangian description, the mesh nodes move exactly with the motion of the fluid or solid and thus are often characterized as moving domain problems. Lagrangian problems are also often described in terms of the reference configuration (the initial computational geometry before any motion) and the current configuration (the current state of the geometry and mesh nodes).ALE combines these descriptions into a formulation that is convenient for describing FSI problems, where the fluid mesh motion does not follow the motion of the fluid, but instead the mesh motion is determined by some moving scheme, which usually includes the solution of equations of elasticity. More information can be found in the literature.

References

[1] Ju Liu, Alison L. Marsden. A unified continuum and variational multiscale formulation for fluids, solids, and fluid–structure interaction. Computer Methods in Applied Mechanics and Engineering; 1 August 2018.

[2] Jie Cheng and Lucy T. Zhang. A General Approach to Derive Stress and Elasticity Tensors for Hyperelastic Isotropic and Anisotropic Biomaterials. International Journal of Computational Methods; Vol. 15, No. 04, 1850028 (2018).

[3] T. Christian Gasser, Ray W Ogden and Gerhard A Holzapfel.Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. Royal Society. 28 September 2005.

[4] J. M. Guccione, A. D. McCulloch, L. K. Waldman. Passive Material Properties of Intact Ventricular Myocardium Determined From a Cylindrical Model.J Biomech Eng. Feb 1991, 113(1): 42-55.

[5] Gerhard A. Holzapfel and Ray W. Ogden. Constitutive modelling of passive myocardium: a structurally based framework for material characterization. Royal Society. 13 September 2009.

[6] Lei Shi, Ian Y. Chen, Hiroo Takayama, Vijay Vedula. An optimization framework to personalize passive cardiac mechanics.Computer Methods in Applied Mechanics and Engineering. 1 December 2024.