A sharp interface method for an immersed viscoelastic solid

TL;DR

This paper introduces a sharp interface method for immersed viscoelastic solids that accurately captures pressure discontinuities, improving convergence and reducing errors in fluid-structure interaction simulations.

Contribution

It presents a novel pressure decomposition approach enabling high-accuracy resolution of pressure discontinuities in immersed boundary-finite element methods for large deformations.

Findings

Enhanced convergence rates for displacements, velocities, stresses, and pressures.

Smaller errors at comparable degrees of freedom.

Validated on analytic solutions, benchmark problems, and active contracting structures.

Abstract

The immersed boundary-finite element method (IBFE) is an approach to describing the dynamics of an elastic structure immersed in an incompressible viscous fluid. In this formulation, there are discontinuities in the pressure and viscous stress at fluid-structure interfaces. The standard immersed boundary approach, which connects the Lagrangian and Eulerian variables via integral transforms with regularized Dirac delta function kernels, smooths out these discontinuities, which generally leads to low order accuracy. This paper describes an approach to accurately resolve pressure discontinuities for these types of formulations, in which the solid may undergo large deformations. Our strategy is to decompose the physical pressure field into a sum of two pressure-like fields, one defined on the entire computational domain, which includes both the fluid and solid subregions, and one defined only on the solid subregion. Each of these fields is continuous on its domain of definition, which enables high accuracy via standard discretization methods without sacrificing sharp resolution of the pressure discontinuity. Numerical tests demonstrate that this method improves rates of convergence for displacements, velocities, stresses, and pressures, as compared to the conventional IBFE method. Further, it produces much smaller errors at reasonable numbers of degrees of freedom. The performance of this method is tested on several cases with analytic solutions, a nontrivial benchmark problem of incompressible solid mechanics, and an example involving a thick, actively contracting torus.

Equations

$$\sigma (\mathbf{x},t)=-\mathbb{I}p(\mathbf{x},t)+\mu (\nabla \mathbf{u}(\mathbf{x},t)+\nabla \mathbf{u}(\mathbf{x},t{)}^T)+\{\begin{array}{ll}{\mathbf{\sigma }}^{\text{e}}(\mathbf{x},t),& \mathbf{x}\in {\Omega }_t^{\text{s}},\\ 0,& \mathbf{x}\in {\Omega }_t^{\text{f}},\end{array}$$

$${\mathbb{P}}^{\text{e}}(\mathbf{X},t)=J(\mathbf{X},t){\mathbf{\sigma }}^{\text{e}}(\mathbf{\chi }(\mathbf{X},t),t){\mathbb{F}}^{-T}(\mathbf{X},t).$$

$$\rho \left(\frac{\partial \mathbf{u}}{\partial t}(\mathbf{x},t)+\mathbf{u}(\mathbf{x},t)\cdot \nabla \mathbf{u}(\mathbf{x},t)\right)$$

$$\nabla \cdot \mathbf{u}(\mathbf{x},t)$$

$$\mathbf{f}(\mathbf{x},t)$$

$$-{\int }_{\partial U}{\mathbb{P}}^{\text{e}}(\mathbf{X},t)\mathbf{N}(\mathbf{X})\delta (\mathbf{x}-\mathbf{\chi }(\mathbf{X},t))dA,$$

$$\frac{\partial \mathbf{\chi }}{\partial t}(\mathbf{X},t)$$

$${\mathbb{P}}^{\text{e}}=\frac{\partial W}{\partial \mathbb{F}}.$$

$$[g(\mathbf{x})]=\underset{\epsilon \to 0}{\lim }g(\mathbf{x}+\epsilon \mathbf{n})-\underset{\epsilon \to 0}{\lim }g(\mathbf{x}-\epsilon \mathbf{n}),\mathbf{x}\in {\Gamma }_t^{\text{fs}}.$$

$$[(\nabla \mathbf{u})\mathbf{t}]=[(\nabla \mathbf{u})\mathbf{b}]=0.$$

$$\frac{d}{ds}{u}_i(\mathbf{\beta }(s))=\frac{d\mathbf{\beta }}{ds}\cdot \nabla u_i(\mathbf{\beta }(s))=\mathbf{t}\cdot \nabla u_i(\mathbf{\beta }(s)).$$

$$\displaystyle\hat{{\bm{x}}}=\begin{pmatrix}n_{x}&n_{y}&n_{z}\\ t_{x}&t_{y}&t_{z}\\ b_{x}&b_{y}&b_{z}\end{pmatrix}({{\bm{x}}}-{\color[rgb]{0,0,0}{{\bm{x}}}^{\star}}):={{\bm{T}}}\,({{\bm{x}}}-{\color[rgb]{0,0,0}{{\bm{x}}}^{\star}}),$$

$$0=[(\nabla \mathbf{u})\mathbf{t}]=[({\nabla }_{\hat{\mathbf{x}}}\mathbf{u})\mathbf{T}\mathbf{t}].$$

$$\left[\frac{\partial u_1}{\partial \hat{y}}\right]=\left[\frac{\partial u_2}{\partial \hat{y}}\right]=\left[\frac{\partial u_3}{\partial \hat{y}}\right]=0.$$

$$\left[\frac{\partial u_1}{\partial \hat{z}}\right]=\left[\frac{\partial u_2}{\partial \hat{z}}\right]=\left[\frac{\partial u_3}{\partial \hat{z}}\right]=0.$$

$$0=[\nabla \cdot \mathbf{u}]=\left[\frac{\partial u_1}{\partial \hat{x}}\right]\frac{\partial \hat{x}}{\partial x}+\left[\frac{\partial u_2}{\partial \hat{x}}\right]\frac{\partial \hat{x}}{\partial y}+\left[\frac{\partial u_3}{\partial \hat{x}}\right]\frac{\partial \hat{x}}{\partial z}.$$

$$\mathbf{n}\cdot [\nabla \mathbf{u}]\mathbf{n}=\mathbf{n}\cdot [{\nabla }_{\hat{\mathbf{x}}}\mathbf{u}]\mathbf{T}\mathbf{n}=\left[\frac{\partial u_1}{\partial \hat{x}}\right]\frac{\partial \hat{x}}{\partial x}+\left[\frac{\partial u_2}{\partial \hat{x}}\right]\frac{\partial \hat{x}}{\partial y}+\left[\frac{\partial u_3}{\partial \hat{x}}\right]\frac{\partial \hat{x}}{\partial z}=0.$$

$$[p]=-\mathbf{n}\cdot {\mathbf{\sigma }}^{\text{e}}\mathbf{n}.$$

$$-[p]\mathbf{n}+\mu ([\nabla \mathbf{u}]+[\nabla {\mathbf{u}}^T])\mathbf{n}-{\mathbf{\sigma }}^{\text{e}}\mathbf{n}=0.$$

$$p(\mathbf{x},t)=\pi (\mathbf{x},t)+\{\begin{array}{ll}\phi (\mathbf{x},t),& \text{if }\mathbf{x}\in {\Omega }_t^{\text{s}},\\ 0,& \text{otherwise.}\end{array}$$

$${\tilde{\mathbb{P}}}^{\text{e}}={\mathbb{P}}^{\text{e}}-J\phi {\mathbb{F}}^{-T}.$$

$${\tilde{\mathbf{\sigma }}}^{\text{e}}=J^{-1}{\tilde{\mathbb{P}}}^{\text{e}}{\mathbb{F}}^T={\mathbf{\sigma }}^{\text{e}}-\phi \mathbb{I}.$$

$$\phi (\mathbf{x},t)=\mathbf{n}\cdot {\mathbf{\sigma }}^{\text{e}}\mathbf{n},\mathbf{x}\in {\Gamma }_t^{\text{fs}}.$$

$$\mathbf{n}da=JdA{\mathbb{F}}^{-T}\mathbf{N}\text{and}{\mathbf{\sigma }}^{\text{e}}\mathbf{n}da={\mathbb{P}}^{\text{e}}\mathbf{N}dA.$$

$$\mathbf{n}=\frac{{\mathbb{F}}^{-T}\mathbf{N}}{\Vert {\mathbb{F}}^{-T}\mathbf{N}\Vert }\text{and}{\mathbf{\sigma }}^{\text{e}}\mathbf{n}=J^{-1}\frac{{\mathbb{P}}^{\text{e}}\mathbf{N}}{\Vert {\mathbb{F}}^{-T}\mathbf{N}\Vert },$$

$$\phi (\mathbf{X},t)=J^{-1}\frac{{\mathbb{F}}^{-T}\mathbf{N}}{\Vert {\mathbb{F}}^{-T}\mathbf{N}{\Vert }^2}\cdot {\mathbb{P}}^{\text{e}}\mathbf{N},\mathbf{X}\in \partial U.$$

$$-{\nabla }_{\mathbf{X}}^2\phi (\mathbf{X},t)$$

$$\phi (\mathbf{X},t)$$

$$\frac{\partial }{\partial t}\phi (\mathbf{X},t)-\gamma {\nabla }_{\mathbf{X}}^2\phi (\mathbf{X},t)$$

$$\phi (\mathbf{X},t)$$

$$\phi (\mathbf{X},0)$$

$$\rho \left(\frac{\partial \mathbf{u}}{\partial t}(\mathbf{x},t)+\mathbf{u}(\mathbf{x},t)\cdot \nabla \mathbf{u}(\mathbf{x},t)\right)$$

$$\nabla \cdot \mathbf{u}(\mathbf{x},t)$$

$$\tilde{\mathbf{f}}(\mathbf{x},t)$$

$$-{\int }_{\partial U}{\tilde{\mathbb{P}}}^{\text{e}}(\mathbf{X},t)\mathbf{N}(\mathbf{X})\delta (\mathbf{x}-\mathbf{\chi }(\mathbf{X},t))dA,$$

Full text

1 Introduction

The immersed boundary (IB) method is a general approach to modeling fluid–structure interaction that was introduced by Peskin [10, 11] for heart valve dynamics. Its strength lies in the representation of the fluid in Eulerian form, which enables approximation on a fixed, Cartesian mesh, and the representation of the solid in Lagrangian form. This approach is appealing because it does not require discretizations that conform to the fluid–structure interface, but instead relies on integral transforms with Dirac delta function kernels to describe interactions between the Eulerian and Lagrangian frames. In practice, conventional IB methods use discretizations of these transforms with regularized delta function kernels. This has the effect of smoothing discontinuities in the pressure and viscous stress that generically appear at the interface between the fluid and solid.

Many techniques have been developed to improve the accuracy in numerical treatments of discontinuities arising in elliptic problems and fluid–structure interaction formulations. One set of approaches focuses on immersed interfaces, i.e., structures that have codimension 1 with respect to the ambient space. LeVeque and Li considered elliptic interface problems and altered finite difference stencils to properly incorporate jump conditions in the discretized equations [6]. We also note the work of Bedrossian et al. in which interface conditions are enforced using a Lagrange multiplier approach [1].

With respect to codimension 1 interfaces immersed in an incompressible viscous fluid, relevant work includes that of Ye et al. [19] for stationary structures and Udaykumar et al. [16], for moving structures. In these methods, the equations of motion are approximated using finite volume methods, and modifications to the discrete operators are developed that account for jump conditions at the immersed interface. Related work by Seo et al. [14] uses a cut–cell approach to obtain higher volume accuracy. Local finite–volume approximations are used in cells that are cut by the interface to modify the equations of motion to account for the imposed jump or boundary conditions. Lee and LeVeque also considered this problem, and imposed jump conditions directly in the Poisson solve for the pressure field used in the discrete solution to the incompressible Navier–Stokes equations [5]. The book by Li and Ito contains a review of some numerical approaches for these types of problems [7].

The focus of this paper is on a sharp interface method for solids with codimension 0. Recent work in high order methods for stationary solids with codimension 0 includes the immersed boundary smooth extension method (IBSE) of Stein et al. [15]. Their approach involves the construction of smooth extensions of solutions from the fluid domain to an extension domain that typically overlaps with the solid region. In the IBSE method, these extensions are spectrally computed from a sequence of harmonic problems and are used to modify forcing functions in the problem formulation. Our approach is similar to the IBSE method in that we solve a harmonic problem and modify the forcing function so the pressure field is continuous. Unlike the approach of Stein et al., however, our approach is able to handle solids undergoing large deformations. We also employ standard finite–volume methods for the Eulerian equations and finite–element methods for the Lagrangian equations.

This study uses the immersed boundary–finite element method (IBFE), in which solid displacements and forces are approximated using finite element discretizations [3]. The first version of such a method appeared in the work of Wang and Liu [18], and was expanded upon in the work of Zhang et al. [20] and Liu et al. [8]. Our paper uses a version of the immersed boundary finite element method based on one introduced by Boffi et al. [2]. Boffi et al. use finite element discretizations for the Eulerian and Lagrangian variables, treat the delta function kernel in a variational way, and obtain suboptimal convergence rates for test problems with respect to the approximation spaces used. In our work, we use the model described by Boffi et al. [2], a finite element approximation for the discretization of the solid, and a second-order staggered–grid finite volume scheme for the discretization of the fluid. This approach was detailed by Griffith and Luo [3].

The method described herein follows the method of Griffith and Luo, but avoids regularizing pressure discontinuities at fluid–structure interface by splitting the physical pressure into the sum of two pressure–like fields that are continuous on their domains of definition. Specifically, we use a global pressure–like field that is defined over the entire Eulerian computational domain along with a second field that is restricted to the Lagrangian domain. Boundary conditions for the Lagrangian pressure–like field ensure that the total pressure recovers the correct jump conditions. This allows us to obtain higher–order accuracy while using only standard discretization methods. Our approach requires one additional harmonic or diffusion equation solve per timestep, compared to the conventional formulation of the IBFE method.

2 Equations of motion and jump conditions

We consider a coupled fluid–structure interaction problem in a domain $\Omega=\Omega_{t}^{\text{s}}\cup\Omega_{t}^{\text{f}}$, with the solid and fluid domains at time $t$ defined to be $\Omega_{t}^{\text{s}}$ and $\Omega_{t}^{\text{f}}$ respectively. The reference configuration of the solid is denoted by $U$. The motion map for the solid is defined to be $\bm{\chi}(\cdot,t)$ so that ${{\bm{x}}}=\bm{\chi}({{\bm{X}}},t)$ is the current position at time $t$ of the reference coordinate ${{\bm{X}}}\in U$. The image $\bm{\chi}(U,t)$ is the current configuration of the solid at time $t$, so that $\Omega_{t}^{\text{s}}=\bm{\chi}(U,t)$ and $\Omega_{t}^{\text{f}}=\Omega\sim\bm{\chi}(U,t)$.

The fluid–structure system that we consider in this paper is described by the total Cauchy stress tensor,

[TABLE]

in which ${\bm{\sigma}}^{\text{e}}$ is the elastic Cauchy stress tensor of the solid, $\mu$ is the dynamic viscosity, $p$ is the pressure, and ${{\bm{u}}}$ is the velocity. Denote the deformation gradient $\mathbb{F}=\partial{{\bm{\chi}}}/{\partial{{\bm{X}}}}$ and its determinant $J=\text{det}(\mathbb{F})$. The first Piola–Kirchoff elastic stress corresponding to ${{\bm{\sigma}}}^{\text{e}}$ is

[TABLE]

As shown by Boffi et al. [2], the equations of motion for this fluid–structure system are:

[TABLE]

in which $\rho$ is the density of the fluid. Equation (1) expresses balance of momentum for the fluid and solid. Incompressibility is imposed in equation (2). Equation (3) explicity defines the force density exerted from the solid onto the fluid in terms of the first Piola–Kirchoff stress. This term includes a volumetric force density and a surface force density, the latter of which generates discontinuities in the pressure, derivatives of the pressure, and derivatives of the velocity at the fluid–structure interface. Equation (4) requires the velocity of the solid to equal the velocity of the background fluid, the so called no–penetration and no–slip condition.

In this paper, we assume the solid is hyperelastic. This implies that the stress is determined by a strain energy density $W$ via

[TABLE]

This assumption is not a requirement of the IBFE approach or our sharp interface method.

We recall some results for deriving a jump condition in the pressure. Sketches of proofs can be found in the work of Lai and Li [4] and Peskin and Printz [12], and we recall them here in our notation for completeness. Let ${{\bm{n}}}={{\bm{n}}}({{\bm{x}}},t)$ and ${{\bm{N}}}={{\bm{N}}}({{\bm{X}}})$ denote the outward unit normal vectors to the solid region in the current and reference configurations respectively. The following notation is useful in discussing jumps of variables at the fluid–structure interface in the current configuration, defined as $\Gamma_{t}^{\text{fs}}=\partial\Omega_{t}^{\text{s}}$.

Definition 1**.**

The jump of a scalar valued function across the fluid–structure interface is defined as:

[TABLE]

The jump of vector or tensor valued variables is defined in the same way, componentwise.

When we use the square bracket notation for the jump described above, we implicitly assume this jump is evaluated at the fluid–structure interface $\Gamma_{t}^{\text{fs}}$. The following results are necessary in deriving a jump condition for the pressure under the assumption the velocity field is continuous.

Lemma 1**.**

Let ${{\bm{t}}}$ and ${{\bm{b}}}$ be the unit tangent vectors at the fluid–structure interface. The tangential derivatives of the velocity are continuous, i.e.

[TABLE]

Proof.

Consider a parametrized curve ${{\bm{\beta}}}={{\bm{\beta}}}(s)$ defined on $\Gamma_{t}^{\text{fs}}$ which contains the point at which we consider the jump. This curve is constructed so its tangent vector $d{{\bm{\beta}}}/ds$ is equal to ${{\bm{t}}}$ at this point. Because the velocity field is continuous, we can consider the velocity evaluated along this curve. The derivative of a component of the velocity field $u_{i}$ along this curve is the tangential derivative,

[TABLE]

This calculation shows the tangential derivative is defined along the boundary and must be continuous. A similar argument is applied for the derivative in the direction of ${{\bm{b}}}$ and for the other components of the velocity field. ∎

Lemma 2**.**

$[{{\bm{n}}}\cdot(\nabla{{\bm{u}}})\,{{\bm{n}}}]=[{{\bm{n}}}\cdot(\nabla{{\bm{u}}})^{T}\,{{\bm{n}}}]=0.$**

Proof.

We follow the argument in [4]. To establish notation, let ${{\bm{x}}}^{\star}$ be the point in the current configuration on $\Gamma_{t}^{\text{fs}}$ where we consider the jump, let ${{\bm{n}}}=(n_{x},n_{y},n_{z})$ be the unit normal vector at ${{\bm{x}}}^{\star}$, and let the orthonormal pair ${{\bm{t}}}=(t_{x},t_{y},t_{z})$ and ${{\bm{b}}}=(b_{x},b_{y},b_{z})$ span the tangent space. Given an arbitrary point ${{\bm{x}}}$, consider a linear transformation to a new point $\hat{{\bm{x}}}$ defined by

[TABLE]

in which $\hat{{\bm{x}}}$ are local coordinates for ${{\bm{x}}}$ expressed in the basis $\{{{\bm{n}}},{{\bm{t}}},{{\bm{b}}}\}$, and ${{\bm{T}}}$ is the matrix containing these vectors in its rows. The notation $\nabla$ refers to derivatives with respect to the physical Cartesian coordinates. At this point, we invoke continuity of tangential derivatives of $\nabla{{\bm{u}}}$ from Lemma 1. By the chain rule, on $\Gamma_{t}^{\text{fs}}$ we have

[TABLE]

With ${{\bm{u}}}=(u_{1},u_{2},u_{3})$ and $\hat{{\bm{x}}}=(\hat{x},\hat{y},\hat{z})$, by orthogonality of the first and last rows of ${{\bm{T}}}$ with respect to ${{\bm{t}}}$, the above statement reads

[TABLE]

The same argument with ${{\bm{b}}}$ gives

[TABLE]

The incompressibility condition, written in terms of $\hat{{\bm{x}}}$, along with (5)–(6) implies

[TABLE]

Finally, we examine

[TABLE]

∎

The discontinuity in the pressure field is derived in the following lemma.

Lemma 3**.**

The pressure satisfies the following jump condition at the fluid–structure interface:

[TABLE]

Proof.

Continuity of the traction vector $[{{\bm{\sigma}}}\,{{\bm{n}}}]=0$ on $\Gamma_{t}^{\text{fs}}$ implies:

[TABLE]

Taking the inner product of this equation with ${{\bm{n}}}$ and using Lemma 2, we obtain the jump condition for the pressure. ∎

3 Description of the method

The method introduced here splits the physical pressure field $p$ into a sum of two components, $\varphi$ and $\pi$. The field $\varphi$ is defined on the solid domain $\Omega_{t}^{\text{s}}$ as a solution to a harmonic problem, to be specified below. The field $\pi$ is defined on $\Omega$ and, by choosing appropriate boundary conditions for $\varphi$, is continuous at the fluid–structure interface. Once $\varphi$ and $\pi$ are known, the physical pressure is recovered by $p=\pi+\varphi$. Explicitly, the physical pressure is defined as:

[TABLE]

This splitting is achieved by modifying the elastic stress ${{\bm{\sigma}}}^{\text{e}}$ so the normal component of the modified traction vanishes on $\Gamma_{t}^{\text{fs}}$. We take the field $\varphi$ to modify the first Piola–Kirchoff stress in the following way:

[TABLE]

Note that $\varphi$ has units of pressure, and we define a modified elastic Cauchy stress $\tilde{{\bm{\sigma}}}^{\text{e}}$ as:

[TABLE]

To encode the pressure discontinuity in the boundary condition for $\varphi$, we require that it satisfy the interface condition:

[TABLE]

This condition removes the normal component of the modified traction at the fluid–structure interface, rendering the pressure–like field continuous; refer to Lemma 3. The pressure–like field $\pi$ is defined to be the pressure solution to the equations of motion when solved with this modified stress. Notice that $\pi$ remains a Lagrange multiplier for $\nabla\cdot{{\bm{u}}}=0$, and its value is completely determined (at least up to an additive constant, depending on boundary conditions) by imposing this constraint.

It is convenient to compute $\varphi$ in the reference configuration of the solid. To express the interface condition for $\varphi$ in the reference configuration, we use Nanson’s relation and the definition of the first Piola–Kirchoff stress,

[TABLE]

These equations imply

[TABLE]

which results in an interface condition for $\varphi$ formulated in the reference configuration:

[TABLE]

We consider two approaches for computing $\varphi$. The first approach is to compute $\varphi$ as a solution to a steady state harmonic problem at each timestep:

[TABLE]

The second approach requires $\varphi$ to solve a diffusion equation each timestep, for some diffusion constant $\gamma$:

[TABLE]

The initial condition $\varphi_{0}$ is taken to be the solution to the harmonic problem at time $t=0$. The parameter $\gamma$ can help to control the condition number for the linear system arising from spatial discretization. We remark that the solutions to either (7)–(8) or (9)–(11) are generally different, resulting in a different splitting of the physical pressure $p=\pi+\varphi$ on the current configuration of the solid.

The equations of motion for this sharp interface formulation of the IBFE scheme, using the modified stress and the steady state harmonic equation for $\varphi$, can be stated as:

[TABLE]

The equations are similar for the case when $\varphi$ is computed via the diffusion equation and are omitted for brevity.

4 Numerical approximation

The numerical approximation for the equations of motion follows [3]. The time step is denoted $\Delta t$. The Cartesian grid spacing parameter is denoted $h$. The mesh factor $M_{\text{fac}}$, referenced in some of the results below, corresponds to the ratio between the approximate edge length in the solid finite element mesh and the Cartesian grid spacing $h$. Solid displacements and forces are approximated using finite elements via the “unified weak formulation” as described by Griffith and Luo [3], for which we seek an approximate volumetric force density ${{\bm{G}}}={{\bm{G}}}({{\bm{X}}},t)$ that is variationally equivalent to the sum of the volumetric ($\nabla_{{\bm{X}}}\cdot\mathbb{P}^{\text{e}}$) and surface ($\mathbb{P}^{\text{e}}\,{{\bm{N}}}$) Lagrangian force densities in (3). More explicitly, given some finite element space $\mathcal{F}_{h}$, the approximate volumetric force density satisfies:

[TABLE]

Nodal Lagrange finite elements are used for the finite element space $\mathcal{F}_{h}$ defined on the solid mesh, and in all simulations, we use bilinear ($Q^{1}$) elements.

The equations for $\varphi$, either the steady state harmonic equation (7)–(8) or the diffusion equation (9)–(11), are discretized in space using a standard finite element discretization. Bilinear elements are also used for the finite element space, and Dirichlet boundary conditions are imposed via the penalty method. The diffusion equation is discretized in time using Crank–Nicolson. The resulting linear systems are solved with GMRES, using an incomplete $LU$ preconditioner. For all numerical experiments, relative linear solver tolerances are set to $1\times 10^{-12}$.

The fluid equations are discretized using a second–order accurate staggered grid method, in which the pressure is approximated at the cell center and the velocity components are approximated at the edges (in two spatial dimensions) or faces (in three spatial dimensions) of the Cartesian grid cell. Lagrangian–Eulerian interaction through the delta function kernels is discretized by designing approximations to these operators that are discrete adjoints. Details are provided in prior work [3].

5 Results

To test this method, we consider four examples. The first two involve a thick ring that is either pre–stressed in its reference configuration or inflated to a steady–state final configuration with a fluid source. Both cases involve substantial pressure discontinuities along the fluid–structure interface, which we demonstrate are well resolved by the method introduced herein. Comparisons to analytic solutions reveal the overall errors in velocities, displacements, stresses and pressures are much smaller with our method, as compared to the conventional IBFE formulation, and converge at faster rates. The third example describes compression of a neo–Hookean block of material. This problem is more challenging because it contains discontinuities in the surface forces that generate solid displacements on the the top and bottom of the block. Further, the material model for the compressed block includes a volumetric energy useful for penalizing compressible deformations [17], but this term prominently contributes to discontinuities at the fluid–structure interface that present additional numerical challenges. The fourth example, inspired by the work of McQueen and Peskin [9], applies the method to an actively contracting thick torus. The contraction produces a transient velocity field and also large pressures in the toroidal wall which leads to pressure discontinuities at the fluid–structure interface.

5.1 Thick orthotropic static ring

For our first example, we consider a two–dimensional thick ring introduced in by Boffi et al. [2]. The Lagrangian curvilinear coordinates of the ring are ${{\bm{s}}}=(s_{1},s_{2})\in U=[0,2\pi R]\times[0,w]$, and its initial configuration is defined by the motion map at $t=0$:

[TABLE]

with $R=0.25$ mm and $w=0.0625$ mm. The constitutive model is taken to be

[TABLE]

in which $\mu_{e}$ is the elastic stiffness coefficient. Note that in this model, the Lagrangian coordinates are not the same as the reference coordinates. This model is initially at equilibrium, and based on the definition of the motion map from curvilinear Lagrangian coordinates to reference coordinates, leads to a discontinuous pressure field at the fluid–structure interface. An analytic solution is available for this problem [2], and we use it to compute convergence rates for the pressure and velocity. Under the condition that the mean of the pressure field is zero, and with $p_{0}=\frac{\pi\mu_{e}}{3w}\left(3wR+R^{2}-\frac{(R+w)^{3}}{R}\right)$ and $r=\|{{\bm{x}}}-(0.5,0.5)\|$, the exact pressure field is:

[TABLE]

For the numerical simulations, the fluid domain is taken to be $[0,L]^{2}$ with $L=1$ mm. In this test we set $\mu_{e}=1$ N/mm, $\rho=$1 kg/mm3 and $\mu=$1 N$\cdot$s/mm2. The velocity is set to zero at the boundary of the fluid domain, and the pressure field $p$ and pressure–like fields $\varphi$ and $\pi$ are normalized to have zero mean. We use the spatial discretization parameter $h=L/N$ with $N=2^{m}$ for $m=5,6,7,8,$ and $9$. The time step size is $\Delta t=0.25\times h$, we set $M_{\text{fac}}=2$, and we perform simulations to a final time of $T=0.01$ s, at which time we compute the errors on the Cartesian grid. Cartesian representations of the pressure–like fields $\pi$ and $\varphi$ computed using the sharp interface method are shown in Figure 1. Notice the continuity of $\pi$ at the fluid–structure interface, because the discontinuity in the physical pressure field $p$ is accounted for in the boundary condition for $\varphi$.

The pressure field computed with the original IBFE method, along with physical pressure field $p=\varphi+\pi$ from the sharp interface method, are shown in Figure 2 on the left and right respectively. The pressure field determined by the original IBFE method contains artifacts at the fluid–structure interface, whereas the sharp interface approach cleanly resolves the pressure discontinuities.

The impact of the sharp interface method in this example can be quantified by considering the velocity and pressure errors. The errors are computed on the Cartesian grid by projecting $\varphi$ onto the grid. This representation of $\varphi$ is visualized on the right in Figure 1. We remark that the errors displayed are absolute, to highlight the differences in absolute error between the original IBFE method and the sharp interface method. Figure 3 shows the errors for the original IBFE method. For velocity, we obtain approximately a rate of 2 for the $L^{1}$ error, a rate of 1.5 for the $L^{2}$ error, and a rate of 1 for the $L^{\infty}$ error. For the pressure, we get a rate of 1 for the $L^{1}$ error, a rate of 0.5 for the $L^{2}$ error, and non–convergence for the $L^{\infty}$ error. This is consistent with the results reported in [3].

Figure 4–6 show the corresponding errors for the different formulations of the sharp interface method. Improved convergence rates are clearly apparent. Each of these figures corresponds to a different method for computing $\varphi$; Figure 4 uses the steady state formulation and Figures 5 and 6 use the diffusion formulation with $\gamma=1$ and $\gamma=h$ respectively. In all cases, the velocity error converges to zero at a second order rate in each norm. The $L^{1}$ error for the pressure converges with a rate of 2, the $L^{2}$ error has a rate of 1.5, and the $L^{\infty}$ error has a rate of 1. Perhaps most significantly, the sharp interface method yields an improvement in the errors of several orders of magnitude.

Figure 7 plots the average number of linear solver iterations required to compute $\varphi$, as we vary $h$. We consider the steady state formulation for $\varphi$ as well as the diffusion formulation, with diffusion constant $\gamma$ set equal to either $1$ or $h$. Notice that the number of iterations scale more mildly for cases that determine $\varphi$ with the diffusion formulation. If $\gamma=h$, we observe that the number of iterations is essentially independent of $h$.

5.2 Thick inflating ring

For this test, we consider the same ring as in the first example, but with a different constitutive model. In this case, the reference configuration is circular, i.e.

[TABLE]

In particular, we note the reference coordinates are the initial coordinates. The material model is defined with respect to the polar deformation gradient $\mathbb{F}_{\text{polar}}$ and polar first Piola–Kirchoff stress $\mathbb{P}^{\text{e}}_{\text{polar}}$ as follows:

[TABLE]

A volume of fluid $A_{\text{add}}$ is injected gradually in the center of the ring, and the pressure is held at zero outside of the ring. At steady state, the exact solution for the pressure field is

[TABLE]

where $r_{\text{in}}=r(R_{\text{in}})$ and $r_{\text{out}}=r(R_{\text{out}})$ are given as

[TABLE]

A derivation can be found in the appendix, along with an explanation of the notation. We remark that the pressure discontinuities can be made to be very large in this example, by increasing the value of $\mu_{e}$. The pressure in the center of the ring is takes a positive value, but within the solid, the pressure takes large negative values for a large choices of $\mu_{e}$. In this test we set $\mu_{e}=10^{4}$ N/mm2, $\rho=$1 kg/mm3 and $\mu=$1 N$\cdot$s/mm2.

For our simulations, the fluid volume $A_{\text{add}}$ is added to the ring center over the first 0.1 s, and we run the model to a final time of 1 s to reach an approximate steady state. A combination of normal traction and tangential no–slip velocity boundary conditions are applied to the computational domain $\Omega$ to hold the pressure at zero outside of the ring. The computational domain $\Omega$ is $[-L,L]^{2}$ with $L=1$ mm. The initial inner and outer radii of the ring are $R_{\text{in}}=0.25$ mm and $R_{\text{out}}=0.3125$ mm respectively. The spatial discretization parameter is $h=L/N$ with $N=2^{m}$ for $m=5,6,7,8,$ and $9$. The time step size is $\Delta t=0.025\times h$, and we set $M_{\text{fac}}=1$.

Figure 8 shows the nontrivial displacement of the solid finite element from its initial to final configuration, with the original IBFE method on the left and the sharp interface method on the right. In this figure, we use the steady state formulation for $\varphi$. Because the problem is axially symmetric, we show only one quarter of the ring. The final configuration is colored, indicating the elementwise value of $J=\det(\mathbb{F})$. In this case, the Eulerian grid is very coarse and contains $32^{2}$ fluid cells; because $M_{\text{fac}}=1$, the initial finite element mesh indicates the approximate size of the fluid cells. Inflation of the ring results in a deformed configuration in which the ring thickness is close to the size of a single Cartesian grid cell. The deformations in the original IBFE method deviate from incompressible deformations (i.e., corresponding to $J=1$), perhaps not surprisingly because the Cartesian grid is coarse with respect to the thickness of the ring in the final configuration. In contrast, the deformations determined by the sharp interface method are much closer to being incompressible, highlighting its utility in numerically conserving volume, at least for coarse finite element meshes.

Figure 9 shows the pressure–like fields $\pi$ and $\varphi$ on the Cartesian grid. As in the static ring case, $\pi$ is continuous at the fluid–structure interface, and the pressure discontinuity is completely accounted for by the boundary conditions imposed in the $\varphi$ field. Figure 10 shows the physical pressure $p$. Results from the original IBFE method are shown on the left and the physical pressure $p=\pi+\varphi$ determined by the sharp interface method are shown on the right. The original IBFE method clearly smears out the pressure discontinuity, whereas the discontinuity is cleanly resolved in the sharp interface method.

Figures 11 and 12 examine the absolute errors in the pressure, velocity, displacement, and elastic stress. The pressure at the center of the ring is sampled using a cosine kernel with radius equal to 0.1 mm. Figure 11 reports results for the original IBFE method. Results for the sharp interface method are shown in Figure 12, with the steady state formulation used compute $\varphi$.

In Figure 11, for the original IBFE method, we observe the same convergence rates for the velocity and pressure as in the static ring test. The velocity errors converge with a rate of 2 in $L^{1}$, a rate of 1.5 in $L^{2}$, and a rate of 1 in $L^{\infty}$, and the pressure errors converge with a rate of 1 in $L^{1}$ and a rate of 0.5 in $L^{2}$. With the conventional method, there is no pointwise convergence of the overall pressure field, and the pressure sampled at the center of the ring converges with a rate of 1. The displacement errors converge with the same rates as the corresponding velocity errors. The elastic stress errors are effectively the errors in the gradient of the displacement, which yields no pointwise convergence of the elastic stress error.

Figures 12 indicates that the sharp interface method gives improved rates of convergence as well as much smaller absolute errors. All velocity and displacement errors converge with a rate of 2, and correspondingly the elastic stress errors converge with a rate of 1. The error rates for the overall pressure field are improved: the $L^{1}$ error converges with a rate of 2, the $L^{2}$ error converges with a rate of 1.5, and the $L^{\infty}$ error converges with a rate of 1. The pressure sampled at the center of the ring converges with a rate of 2. The results for the case where $\varphi$ is computed with the diffusion formulation with $\gamma=1$ are the same. When $\gamma=h$, the results are essentially the same except for the velocity errors. In this case, the velocity errors appear to converge with a rate of 1 in all norms.

5.3 Compressed block

Our third example is related to a standard solid mechanics test problem: compression of a two–dimensional block of isotropic material [13]. The strain energy density is neo–Hookean and is split into isochoric $W^{\text{iso}}$ and dilational $W^{\text{dil}}$ parts. This decomposition is achieved by using the first invariant $I_{1}=\text{tr}(\mathbb{F}^{T}\mathbb{F})$, which is modified by scaling the deformation gradient so that its determinant is 1:

[TABLE]

Then the strain energy density is defined by

[TABLE]

for shear modulus $\mu_{e}$ and bulk modulus $\lambda=\lambda(\nu)$, taken to be a function of a numerical Poisson ratio $\nu$ which is detailed in [17]. The dilational energy $W^{\text{dil}}$ vanishes at the continuous level because $J=1$. In the discretized equations, this term is not necessarily equal to zero, and it can be interpreted as a stabilization that enforces discrete incompressibility [17]. The first Piola–Kirchoff stress takes the form:

[TABLE]

with $\lambda$ related to $\nu$ by:

[TABLE]

As suggested in Ref. [17], we consider values for the numerical Poisson ratio $\nu=-1,0,0.4$.

In this test, we include surface force densities on the boundary of the block to weakly impose Dirichlet boundary conditions for the displacement. First, define the following surface force density which tethers points to their reference configuration:

[TABLE]

Let ${{\bm{e}}}_{1}$ and ${{\bm{e}}}_{2}$ denote the standard Cartesian unit vectors. The reference configuration of the block is setup so its sides are parallel to the Cartesian axes. The surface force density on the top is then defined as

[TABLE]

and on the bottom is defined as

[TABLE]

In particular, the surface force density on the bottom of the block allows the bottom boundary to slide in the horizontal direction. The parameter used to determine the strength of the force is:

[TABLE]

For this example, we consider both discontinuous and continuous loading pressures to compress the block. A discontinuous loading pressure was part of the original benchmark presented in Ref. [13]; our sharp interface method improves results in this case, but it also encounters some difficulties. In this light, we also consider a smoothed loading pressure that is defined to be the discontinuous pressure multiplied by a mollifier.

5.3.1 Discontinuous loading pressure

Define ${{\bm{X}}}=(X_{1},X_{2})$. The discontinuous loading pressure takes the form

[TABLE]

with

[TABLE]

In the reference configuration, the block geometry spans the horizontal direction from $X_{1}=0$ mm to $X_{1}=20$ mm and the vertical direction from $X_{2}=0$ mm to $X_{2}=10$ mm. We set $\tilde{a}=5$ mm, $\tilde{b}=15$ mm, and the loading pressure is applied to the top of the block at $X_{2}=10$ mm. The maximum loading pressure is taken to be $P_{\text{max}}=200$ N/mm2, and we set $\mu_{e}=80.194$ N/mm2. The fluid density is $\rho=$1 kg/mm3 and the viscosity is $\mu=$0.16 N$\cdot$s/mm2.

The loading pressure is ramped up over the first 10 s of the simulation, so $t_{\text{load}}=$ 10 s. The model is run to a final time of 50 s to reach an approximate steady state. As in the inflating ring test, a combination of normal traction and tangential no–slip velocity boundary conditions are applied to the computational domain $\Omega$ to hold the pressure at zero. The computational domain $\Omega$ is $[0,L]^{2}$ with $L=30$ mm, and the block is positioned in its center. The spatial discretization parameter is $h=L/N$ with $N=2^{m}$ for $m=4,5,6,$ and $7$. The time step size is $\Delta t=0.005\times h$ and we set $M_{\text{fac}}=1$.

Figure 13 shows the displaced mesh for the compressed block with the discontinuous loading pressure. For this result, the steady state formulation for $\varphi$ is used and $\nu=0$. Figure 14 examines slices of the pressure field down the center of the block at $X_{2}=15$ mm, with the original IBFE method on the top row and the sharp interface method on the bottom row. As $N$ increases, the pressure slices converge to a sharp profile in both methods. The sharp interface method converges faster and more cleanly resolves pressure discontinuities on the top and bottom faces of the block. For the case $\nu=0.4$, the discontinuities in the pressure field are less clear with the sharp interface method. This observation indicates some subtle interplay between the volumetric energy and the splitting of the pressure field, which requires additional investigation. We consider this issue in the next section with the continuous loading pressure.

In Figure 15, we plot $\pi$ and $\varphi$ for $\nu=0$. The $\pi$ field appears to contain discontinuities at the approximate location of the jumps in the pressure loading force. The sharp interface method presented here assumes that $\pi$ is continuous, so this example presents some challenges for our approach. Even so, the physical pressure field from the sharp interface method, in Figure 16 on the right, more sharply captures pressure discontinuities compared to the original IBFE method.

5.3.2 Continuous loading pressure

The smoothed version of the loading pressure we use takes the form:

[TABLE]

The parameters $a=4$ mm, $b=16$ mm, and all other numerical and physical parameters are the same as in the previous section.

Figure 17 depicts the displaced mesh after the continuous loading pressure is applied, compared to the initial mesh in the background. First, we more closely examine the impact of the volumetric energy on the sharp interface method by displaying the $\pi$ and $\varphi$ fields separately in Figures 18 and 19. When the volumetric energy is zero, i.e. $\nu=-1$, the $\pi$ field has a large gradient at the bottom of the block. When the volumetric energy is included in the strain energy density, corresponding to the cases $\nu=0$ and $\nu=0.4$, the gradients in the $\pi$ field are more mild, with larger values toward the top of the block where the pressure load is applied. The $\varphi$ field in all cases takes it largest values at the top of the block, and in the cases with nonzero volumetric energy has smaller discontinuities on the block boundary.

In Figure 20, we examine the $\pi$ and $\varphi$ fields for the case $\nu=0$. Not surprisingly, the largest pressure discontinuities in this test appear on the top of the block where the loading pressure is applied. Further, the computed $\pi$ field is more smooth than in the case with the discontinuous loading pressure. Figure 21 shows the physical pressure field; on the left are the results from the original IBFE method and on the right is the sum $\pi+\varphi$ from the fields in Figure 20. This sharp interface method better resolves the pressure discontuities, especially in the center of the top of the block where the loading pressure is applied.

Figure 22 examines a slice of the pressure field taken vertically at $X_{2}=15$ mm, for different mesh sizes and different values for $\nu$. We only display slices for the sharp interface method with the steady state formulation for $\varphi$ since results for the diffusion formulation look similar.

Figures 23 and 24 display approximate absolute pressure and velocity errors, computed on the finite element mesh. The errors are computed as the difference between coarse and fine solutions at consecutive levels of refinement. For the original IBFE method, Figure 23, the $L^{1}$ error for the velocity appears to be converging with a rate of 1 and the $L^{\infty}$ error is not converging. The $L^{2}$ converges with a rate less than 1. The pressure errors show a similar trend. The approximate errors for the sharp interface method, shown in Figure 24, are improved. All velocity and pressure errors converge with approximately a rate of 1. Results are similar when using the diffusion equation for $\varphi$ with either $\gamma=1$ or $\gamma=h$.

5.4 Actively contracting thick torus

In our final example, we apply the method to a thick, actively contracting torus which creates a nonzero velocity field. This problem is inspired by the work of McQueen and Peskin [9]. The constitutive model contains an active stress component describing a wave of contraction; this wave pushes fluid within the torus, creating a peristaltic pump.

The geometrical parameters are $R_{1}$, $R_{2}$, and $w$, and the thick toroidal geometry is parametrized with respect to the coordinates ${{\bm{s}}}=(s_{1},s_{2},s_{3})=[-\pi R_{1},\pi R_{1})\times[0,w]\times[-\pi R_{2},\pi R_{2})$. The reference coordinates ${{\bm{X}}}=(X_{1},X_{2},X_{3})$ are defined as:

[TABLE]

We identify this parametrization by a function $\mathcal{G}$, i.e. ${{\bm{X}}}=\mathcal{G}({{\bm{s}}})$. The parameter $w=0.125$ mm is the thickness of the wall, built from inner and outer toroidal surfaces. The parameters $R_{1}=0.25$ mm and $R_{2}=1$ mm are the minor and major radii of the inner toroidal surface respectively.

The elastic stress is a sum of two components, a passive stress and an active stress which depends on a fiber vector field ${{\bm{f}}}_{0}$ specified in the reference configuration:

[TABLE]

The size, strength, and location of the contraction is determined by a tension function $T=T({{\bm{s}}},t)$. The support of this function has measure $\ell_{\text{c}}$ in the third component $s_{3}$ and translates with velocity $\nu_{\text{c}}$. Explicitly, the tension function is defined as

[TABLE]

with $a=-\frac{1}{2}\ell_{\text{c}}$ and $b=\frac{1}{2}\ell_{\text{c}}$. Notice the tension function smoothly decays to zero. The maximum of the tension is ramped up linearly for an amount of time $t_{\text{ramp}}$ to a value $T_{\text{max}}$:

[TABLE]

The finite element mesh for the torus, with the fiber field ${{\bm{f}}}_{0}$ superimposed, is shown in Figure 25. The angle of the fiber field rotates from $60^{\text{o}}$ to $-60^{\text{o}}$ from the outer surface to the inner surface. The mesh is centered in the computational domain defined to be $\Omega=[-2L,2L]\times[-L,L]\times[-2L,2L]$ with $L=0.75$. This domain is discretized with $2N$ cells in the $x$ and $z$ directions and $N$ cells in the $y$ direction, with $N=64$. The finite element mesh is defined so $M_{\text{fac}}$ is approximately 2. Boundary conditions for the computational domain are set to hold the pressure at zero. The timestep size is $\Delta t=0.0025\times h$, with $h=L/N$.

The displacement of the torus during contraction is shown in Figure 26, with time increasing from left–to–right, top–to–bottom. The color indicates the value of the tension function. The torus is not pre–stressed, but the active stress term results in large pressure discontinuities at the fluid–structure interface and induces a nonzero velocity field. Figures 27 and 28 show the velocity field on a slice through the center of the displaced torus for the original IBFE method and the sharp interface method respectively. The Reynolds number is approximately 195, and both methods qualitatively produce the same velocity field. Notice the region of forward flow directly in front of the contracted part of the torus.

Figures 29 and 30 show a slice of the pressure field from the original IBFE method and the sharp interface method respectively. The active stress term creates large pressures in the wall of the torus, leading to pressure discontinuities at the fluid–structure interface. These pressure jumps are smoothed out with the original approach. The sharp interface method resolves these discontinuities.

6 Conclusion

This paper describes a numerical method for resolving pressure discontinuities in an immersed boundary finite element formulation. This method applies to immersed solids of codimension 0 which may undergo large deformations. The pressure is split into two pressure–like fields, one of which is continuous and defined on the entire domain, and the other is defined only on the finite element mesh. Calculation of the pressure–like field defined on the finite element mesh is done by solving either a steady state harmonic problem or a diffusion equation, with boundary conditions that depend on the normal component of the elastic traction.

The method is tested on four examples. One example statically contains pressure discontinuities, and the other three involve large deformations. In each test, the method resolves the pressure discontinuities on the fluid–structure interface. In the thick ring examples with analytic solutions, we see improved convergence for the pressure, velocity, displacement, and elastic stress, and a substantial decrease in the magnitude of errors for modest numbers of degrees of freedom. Further, the inflating ring example demonstrates that the method helps to enforce solid incompressibility, at least for coarse finite element meshes. In the more complex compressed block example, the method more sharply resolves the pressure field for both smooth and discontinuous loading forces. Results with the discontinuous loading force suggest that care must be taken in managing discontinuous forcing on the fluid–structure interface. In the case of a smooth loading force, extrapolated errors for both the pressure and velocity appear to be converging. This example highlights the method’s flexibility in dealing with additional surface forces and volumetric energies with modestly large stabilization parameters [17]. In the final example, we demonstrated the method works in three dimensions with transient flow dynamics and a stress function which depends on a fiber vector field. The active stress term creates large pressures in the toroidal wall, leading to pressure discontuities at the fluid–structure interface which are easily resolved with our approach.

As shown in the compressed block example, the proposed method appears to struggle with large volumetric stabilization parameters. A related question is if this approach can serve as a replacement for the volumetric stabilization technique, given its performance in enforcing incompressible solid deformations for the inflating ring. These questions are topics of future work.

Additional future work includes the resolution of other jump conditions, including the normal derivative of the velocity at the fluid–structure interface, possibly by solving higher order harmonic problems on the solid finite element mesh. Further, a modification of the method may work for codimension 1 structures which correspond to the boundary of a codimension 0 structure. A similar technique may be possible if a volumetric mesh for the corresponding codimension 0 structure is available; in this case, a harmonic problem could be solved on the associated volumetric mesh to impose known jump conditions across the thin interface.

7 Acknowledgments

Charles Puelz was supported in part by the Research Training Group in Modeling and Simulation funded by the National Science Foundation via grant RTG/DMS–1646339. Boyce Griffith was supported in part by NSF DMS 1664645, NSF OAC 1652541, NSF OAC 1450327. The authors thank Charles Peskin for many insightful conversations.

Appendix A Exact solution for thick inflating ring

Consider a two–dimensional thick ring with with inner radius $R_{\text{in}}$ and outer radius $R_{\text{out}}$. We inject fluid on the interior of the ring with volume $A_{\text{add}}$. By incompressibility, the deformed radius within the ring is

[TABLE]

This uses a motion map in polar coordinates, which may be expressed as:

[TABLE]

The deformation gradient in polar coordinates is:

[TABLE]

Let us consider the following material model:

[TABLE]

Assuming incompressibility, $J=1$, the elastic Cauchy stress is

[TABLE]

By force balance, at steady state, we have:

[TABLE]

This implies that

[TABLE]

If we hold the pressure on the outside of the ring equal to zero, we can determine this constant from the jump condition on the outer radius. More precisely:

[TABLE]

Choosing $p^{+}=0$ implies that $p^{-}=-\mu_{e}\frac{A_{\text{add}}}{\pi r_{\text{out}}^{2}}$. But since

[TABLE]

we have determined the constant to be $-\frac{\mu_{e}}{2}\frac{A_{\text{add}}}{\pi r_{\text{out}}^{2}}$. Thus, we obtain

[TABLE]

The jump condition for the pressure on the interior radius allows us to derive a formula for the pressure of the fluid interior to the ring. Similarly, we have

[TABLE]

Solving for $p^{+}$ gives us a formula for the pressure of the fluid interior to the ring. The pressure is:

[TABLE]

As a remark, the implementation of this model uses Cartesian coordinates, and care is needed when converting between Cartesian and polar coordinate systems. Define the reference polar coordinates to be ${{\bm{P}}}=(R,\Theta)$ and the current polar coordinates to be ${{\bm{p}}}=(r,\theta)$. Converting the polar Cauchy stress ${{\bm{\sigma}}}_{\text{polar}}^{\text{e}}$ to the Cartesian version ${{\bm{\sigma}}}_{\text{cart}}^{\text{e}}$, with $x=r\cos\theta$ and $y=r\sin\theta$, uses the Givens rotation:

[TABLE]

Also, the chain rule can be used to convert the Cartesian deformation gradient $\mathbb{F}_{\text{cart}}$ to the polar deformation gradient $\mathbb{F}_{\text{polar}}$:

[TABLE]