A conservative and non-dissipative Eulerian formulation for the simulation of soft solids in fluids

TL;DR

This paper introduces a conservative, non-dissipative Eulerian finite volume method using the reference map technique for simulating incompressible hyperelastic solids in fluids, improving robustness and interface smoothness without artificial viscosity.

Contribution

The authors develop a robust Eulerian formulation with modifications that eliminate artificial viscosity and interface wrinkles, enhancing simulation accuracy for soft solids in fluids.

Findings

Elimination of artificial viscosity in solid regions.

Reduction of interface wrinkles and improved smoothness.

Enhanced numerical robustness and stability.

Abstract

Soft solids in fluids find wide range of applications in science and engineering, especially in the study of biological tissues and membranes. In this study, an Eulerian finite volume approach has been developed to simulate fully resolved incompressible hyperelastic solids immersed in a fluid. We have adopted the recently developed reference map technique (RMT) by Valkov et. al (J. Appl. Mech., 82, 2015) and assessed multiple improvements for this approach.These modifications maintain the numerical robustness of the solver and allow the simulations without any artificial viscosity in the solid regions (to stabilize the solver). This has also resulted in eliminating the striations ("wrinkles") of the fluid-solid interface that was seen before and hence obviates the need for any additional routines to achieve a smooth interface. An approximate projection method has been used to project the velocity field onto a divergence free field. Cost and accuracy improvements of the modifications on the method have also been discussed.

Tables

Table 1. Table 1: Comparison between the present method and the original RMT method by Valkov et al. ( 2015 ) .

	The original RMT (Valkov et al., 2015)	Present approach
Grid	Staggered	Collocated
Nature of the solver	Compressible	Incompressible
Reference map extrapolation	PDE based	Least-squares based
Discrete momentum conservation (inertial terms)	No	Yes
Discretization stencil	One-sided (artificial damping)	Central
Smoothing routines	Required (artificial damping)	No
Global damping	Needed (artificial damping)	Not needed

Table 2. Table 2: Comparison of the error in the cost effective least-square extrapolation procedure vs the hyperbolic PDE extrapolation procedure of Aslam ( 2004 ) . E ξ → ⋅ i ^ subscript 𝐸 ⋅ → 𝜉 ^ 𝑖 E_{\vec{\xi}\cdot\hat{i}} , E ξ → ⋅ j ^ subscript 𝐸 ⋅ → 𝜉 ^ 𝑗 E_{\vec{\xi}\cdot\hat{j}} represents the L 2 subscript 𝐿 2 L_{2} norm error of the x 𝑥 x and y 𝑦 y components of ξ → → 𝜉 \vec{\xi} field computed after one full rotation of the Zalesak disk on a 100 × 100 100 100 100\times 100 grid.

	Least-squares approach	PDE approach
$E_{\overrightarrow{\xi }.\hat{i}}$	$6.72\times {10}^{-9}$	$8.32\times {10}^{-4}$
$E_{\overrightarrow{\xi }.\hat{j}}$	$7.34\times {10}^{-9}$	$4.67\times {10}^{-4}$

Table 3. Table 3: Comparison of the normalized L 2 subscript 𝐿 2 L_{2} norm error E ξ ​ n ​ o ​ r ​ m subscript 𝐸 𝜉 𝑛 𝑜 𝑟 𝑚 E_{\xi norm} for the case of asymmetric star-looking object for various grid sizes.

Grid size	$E_{\overrightarrow{\xi }.\hat{i}norm}$	$E_{\overrightarrow{\xi }.\hat{j}norm}$
$50\times 50$	$1.28\times {10}^{-12}$	$1.35\times {10}^{-12}$
$100\times 100$	$6.39\times {10}^{-13}$	$6.77\times {10}^{-13}$
$200\times 200$	$3.53\times {10}^{-13}$	$3.70\times {10}^{-13}$
$400\times 400$	$4.52\times {10}^{-13}$	$4.59\times {10}^{-13}$

Equations

$$\vec{\xi }(\vec{x},t)=\vec{X},$$

$$\frac{D\vec{\xi }(\vec{x},t)}{Dt}=0.$$

$$\frac{\partial \vec{\xi }(\vec{x},t)}{\partial t}+\vec{u}.\vec{\nabla }\vec{\xi }(\vec{x},t)=0.$$

$$\mathbb{F}(\vec{X},t)=\partial \vec{x}/\partial \vec{X}=[\overrightarrow{\nabla }\overrightarrow{\xi }(\vec{x},t){]}^{-1},$$

$$\frac{\partial (\rho \vec{u})}{\partial t}+\vec{\nabla }.(\rho \vec{u}\otimes \vec{u})=\vec{\nabla }.\underline{\underline{\sigma }},$$

$$\frac{\partial \rho }{\partial t}+\vec{\nabla }.(\vec{u}\rho )=0,$$

$${\underline{\underline{\sigma }}}^s=(det\mathbb{F}{)}^{-1}\mathbb{F}\frac{\partial \bar{\psi }(\mathbb{E})}{\partial \mathbb{E}}{\mathbb{F}}^T-\lambda \text{mathds}1=2(det\mathbb{F}{)}^{-1}\mathbb{F}\frac{\partial \hat{\psi }(\mathbb{C})}{\partial \mathbb{C}}{\mathbb{F}}^T-\lambda \text{mathds}1,$$

$$\hat{\psi }(\mathbb{C})={\mu }^s(tr\mathbb{C}-3),$$

$$\frac{\partial \hat{\psi }(\mathbb{C})}{\partial \mathbb{C}}={\mu }^s\text{mathds}1.$$

$${\underline{\underline{\sigma }}}^s=2{\mu }^s[(\vec{\nabla }\vec{\xi }{)}^{-1}(\vec{\nabla }\vec{\xi }{)}^{-T}]-P\text{mathds}1=2{\mu }^s[(\vec{\nabla }\vec{\xi }{)}^T(\vec{\nabla }\vec{\xi }){]}^{-1}-P\text{mathds}1.$$

$$\underline{\underline{\sigma}}^{s}=2\mu^{s}\begin{bmatrix}(\frac{\partial\alpha}{\partial y})^{2}+(\frac{\partial\beta}{\partial y})^{2}&-\Big{[}(\frac{\partial\alpha}{\partial x})(\frac{\partial\alpha}{\partial y})+(\frac{\partial\beta}{\partial x})(\frac{\partial\beta}{\partial y})\Big{]}\\ -\Big{[}(\frac{\partial\alpha}{\partial x})(\frac{\partial\alpha}{\partial y})+(\frac{\partial\beta}{\partial x})(\frac{\partial\beta}{\partial y})\Big{]}&(\frac{\partial\alpha}{\partial x})^{2}+(\frac{\partial\beta}{\partial x})^{2}\end{bmatrix}-P\mathds{1},$$

$${\underline{\underline{\sigma }}}^f={\mu }^f\left[\left(\vec{\nabla }\vec{u}\right)+{\left(\vec{\nabla }\vec{u}\right)}^T\right]-P\text{mathds}1=\mu \left[\begin{array}{cc}2\frac{\partial u}{\partial x}& \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\\ & \\ \frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}& 2\frac{\partial v}{\partial y}\end{array}\right]-P\text{mathds}1,$$

$$\frac{\partial \vec{\xi }(\vec{x},t)}{\partial t}+\vec{\nabla }\cdot [\vec{u}\vec{\xi }(\vec{x},t)]=0,$$

$$\frac{{\rho }^{n+1}{\vec{u}}_P^{\ast }-{\rho }^n{\vec{u}}_P^n}{\Delta t}=-\vec{\nabla }.{\left(\rho {\vec{u}}_P{\vec{u}}_P\right)}^n,$$

$$\frac{{\rho }^{n+1}{\vec{u}}_P^{\ast \ast }-{\rho }^{n+1}{\vec{u}}_P^{\ast }}{\Delta t}=\vec{\nabla }.\underline{\underline{\sigma }}({\mu }^s,{\mu }^f,{\vec{u}}_P^{\ast },\vec{\xi }).$$

$${\vec{u}}_f^{\ast \ast }={⟨{\vec{u}}_P^{\ast \ast }⟩}_{P\to f}-\left\{\frac{\Delta t}{{\rho }_f^{n+1}}\left[{\left(\vec{\nabla }P\right)}_f^n-F_f^{n+1}\right]\right\},$$

$$\vec{\nabla }.\left[\frac{{\left(\overrightarrow{\nabla \delta }P\right)}_f^{n+1}}{{\rho }_f^{n+1}}\right]=-\vec{\nabla }.\left(\frac{{\vec{u}}_f^{\ast \ast }}{\Delta t}\right),$$

$$P^{n+1}=P^n+\delta P^{n+1}.$$

$${\vec{u}}_f^{n+1}={\vec{u}}_f^{\ast \ast }-\left[\frac{\Delta t}{{\rho }_f^{n+1}}{\left(\overrightarrow{\nabla \delta }P\right)}_f^{n+1}\right].$$

$${\vec{u}}_P^{n+1}={\vec{u}}_P^{\ast \ast }-\Delta t{⟨\frac{(\vec{\nabla }P{)}_f-F_f}{{\rho }_f}⟩}_{f\to P}^{n+1},$$

$$\frac{\partial \alpha }{\partial x}=\frac{{\alpha }_{i+1,j}-{\alpha }_{i-1,j}}{2\Delta x},$$

$${\vec{u}}_f^{\ast \ast }={⟨{\vec{u}}_P^{\ast \ast }⟩}_{P\to f}-\left\{\frac{\Delta t}{{\rho }_f^{n+1}}\left[{\left(\vec{\nabla }P\right)}_f^n-F_f^{n+1}\right]\right\},$$

$$u_{i+1/2,j}=\frac{u_{i,j}+u_{i+1,j}}{2}$$

$$v_{i,j+1/2}=\frac{v_{i,j}+v_{i,j+1}}{2}$$

$$\varphi (\vec{x},t)=\varphi [\vec{\xi }(\vec{x},t),t=0].$$

$$\frac{\partial \vec{\xi }(\vec{x},t)}{\partial t}+H(\vec{x})\vec{\nabla }\cdot [\vec{u}\vec{\xi }(\vec{x},t)]=0,$$

$$H(\vec{x})=\Bigg{\{}\begin{aligned} 1&\hskip 28.45274pt\Omega_{S}\\ 0&\hskip 28.45274ptelse.\end{aligned}$$

$$flag1=\Bigg{\{}\begin{aligned} &1\hskip 28.45274pt\Omega_{S}\\ &0\hskip 28.45274ptelse\end{aligned}$$

$$flag2=\Bigg{\{}\begin{aligned} &1\hskip 28.45274pt\Omega_{S}+\Omega_{E}\\ &0\hskip 28.45274ptelse\end{aligned}$$

$$\underline{\underline{\sigma }}=\hat{H}[\hat{\varphi }(\vec{x},t)]{\underline{\underline{\sigma }}}^f+\left\{1-\hat{H}[\hat{\varphi }(\vec{x},t)]\right\}{\underline{\underline{\sigma }}}^s,$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Soft solids in fluids are ubiquitous in nature. Study of these systems is of practical relevance in science and engineering, especially in the field of biomedicine (Turitto et al., 1972; Wootton and Ku, 1999; Andrews and Low, 1999; Fogelson and Guy, 2004). Some of such applications involve the study of the interaction between micro-bubble collapse-induced shock waves with the tissue in an animal body (Adami et al., 2016), study of the electroporation phenomenon (Neumann et al., 1982), study of hemodynamics and suspension of blood cells (Pozrikidis, 2003, 2010).

Numerical methods to simulate a fluid-solid coupled system, also known as a fluid-structure interaction (FSI) problem can be broadly classified into mesh-based methods and meshfree methods. Further, the mesh-based methods can be subdivided into (a) fully-Eulerian approach, where both fluids and solids are solved on an Eulerian grid, (b) mixed Lagrangian Eulerian approach, where typically fluids are solved on an Eulerian grid and solids are represented using a Lagrangian grid, (c) fully Lagrangian approach, where both fluids and solids are solved using a Lagrangian grid. A detailed classification of various methods used to study FSI problems is shown in the Figure 1.

FSI has historically been studied using a partitioned-based mixed-Lagrangian-Eulerian approaches, where fluid and solid regions are solved separately on different meshes using different methods (see the, arbitrary Lagrangian-Eulerian (ALE) approach of Hu et al. (2001); Hirt et al. (1974); Nitikitpaiboon and Bathe (1993); Hughes et al. (1981); Belytschko (1980), deforming-spatial-domain/stabilized-space-time approach (DSD/SST) of Tezduyar et al. (1992); Hughes and Stewart (1996)). These methods have been widely used to study problems such as flapping wings (Mittal and Tezduyar, 1995; Takizawa et al., 2012), fluid-particle interaction (FPI) (Mittal and Tezduyar, 1994; Johnson and Tezduyar, 1996, 1997b, 1997a, 1999, 2001), patient-specific arterial modeling of cerebral aneurysms (Torii et al., 2004, 2006a, 2006b, 2007a, 2007b, 2008, 2009, 2010a, 2010b, 2011), parachute modeling (Kalro and Tezduyar, 2000; Stein et al., 2000, 2001) wind-turbine rotor aerodynamics (Takizawa et al., 2011b, a), moving hyperelastic particles (Gao and Hu, 2009) and modeling flow in the heart (Watanabe et al., 2004). However, these methods found success mostly in the stiff limit of the solids (Hu, 1996; Johnson and Tezduyar, 1997a) and was found to be too cumbersome for highly deforming solids, since it requires generating a new grid at each time step.

To overcome the cost of partitioned-based approaches for highly deforming solids, monolithic-based solvers were developed, where a single system of equations are solved simultaneously in a coupled manner (Hübner et al., 2004; Michler et al., 2004; Ryzhakov et al., 2010). To account for the effect of presence of solids, a force term that is computed based on the structural configuration of the solid is added to the fluid equations. These methods have been applied to study problems such as modeling rigid particles (Yuki et al., 2007), modeling flexible bodies (Mori and Peskin, 2008; Zhao et al., 2008), red blood cell (Mori and Peskin, 2008; Eggleton and Popel, 1998; Gong et al., 2009). Further, monolithic-based solvers can be sub-divided into (i) fictitious domain (FD) method / distributed Lagrange multiplier method (Glowinski et al., 1999, 2001; Patankar, 2001), (ii) immersed-boundary methods, where the solid region is represented as a boundary (Mittal and Iaccarino, 2005), (iii) immersed-domain method, where the solid region is represented as a body with a finite volume, for example the immersed finite-element method (Liu et al., 2006, 2007; Zhang and Gay, 2007; Wang and Zhang, 2013) and the immersed continuum method of (Wang, 2006, 2007, 2010)

Immersed-boundary methods are known to be the simplest of all the methods. For example, the classical immersed-boundary (IB) method (Peskin, 1972, 1982a, 2002; Kim and Choi, 2006; Huang and Sung, 2009), the direct-forcing method (Mohd-Yusof, 1997; Luo et al., 2007; Mark and van Wachem, 2008; Guy and Hartenstine, 2010), the penalization approach (Kim and Peskin, 2007), the ghost-cell method, the cut-cell finite volume approach (Clarke et al., 1986) and the immersed-interface method (LeVeque and Li, 1994; Li and Lai, 2001; Li et al., 2003; Li and Ito, 2006; Layton, 2009) all use an Eulerian grid for the fluid region and the boundary of the deforming solid is considered as a forcing term in the fluid equations either in the continuous form (continuous forcing methods) or in discrete form (discrete forcing methods). Owing to its simplicity, these methods have been used to study a wide variety of problems such as magnetohydrodynamics of liquid metals (Grigoriadis et al., 2009), complex flows in irregular domains (Fadlun et al., 2000; Udaykumar et al., 1996; Le et al., 2008; Kim et al., 2001; Iaccarino et al., 2003), turbulent flows (Yang and Balaras, 2006), modeling cochlea (Beyer Jr, 1992), modeling flexible fibers (Stockie and Green, 1998; Wang and Layton, 2009), rigid bodies (Mittal et al., 2004), flow-induced vibration (Mittal et al., 2003), biomimetic flight mechanism (Mittal et al., 2002), flow past an airfoil (Ghias et al., 2004), flexible filaments (Zhu and Peskin, 2003), modeling mechanics of heart (Peskin, 1982b; Griffith, 2005), sperm motility near boundaries (Fauci and McDonald, 1995), microswimmers (Dillon et al., 1995) and ship hydrodynamics (Weymouth et al., 2006; Weymouth, 2008) where it’s been coupled with two-fluid solvers. However, these methods are known to not capture the realistic structural response of the solid and often use a linear theory (infinitesimal strain theory) to approximate the stresses and the deformation of the solid.

On the other hand, relatively less popular class of methods to solve FSI problems are the fully Lagrangian approach such as a particle-finite-element method (PFEM) and meshfree methods such as the reproducing kernel particle method (RKPM) and smoothed-particle hydrodynamics (SPH). These methods have been used for applications such as large number of floating bodies in fluid, bed erosion etc. (Onate et al., 2008).

Finally, the relatively newer class of methods are the fully Eulerian approaches. These methods typically use an interface-capturing method that was initially developed to track material interfaces in two-fluid flows (see Mirjalili et al. (2017)). These approaches are inherently cost effective due to a fixed mesh and results in a easily parallelizable computer programs and is particularly very advantageous compared to other methods for highly deforming solids. Some of previous applications of these methods are modeling linear elastic materials (Xiao, 1999), elasto-plastic materials (Udaykumar et al., 2003; Okazawa et al., 2007) and neo-hookean materials (Liu and Walkington, 2001; van Hoogstraten et al., 1994; Dunne, 2006; Sugiyama et al., 2010), modeling flow in complex domains (Nagano et al., 2010). One of the earliest known Eulerian approach that can solve solid regions using “true nonlinear solid constitutive laws” coupled with fluid flow is the Eulerian Godunov method of Miller and Colella (2001)). This has been applied to study mostly elastic-plastic solids (Barton et al., 2010; Ghaisas et al., 2018). However, the main disadvantage of this method is that the method is limited to unbounded domains.

Other well-known fully Eulerian approaches are the level-set method (Cottet et al., 2008; Cottet and Maitre, 2016) and the volume-of-fluid (VOF) method (Sugiyama et al., 2011; Ii et al., 2011, 2012). These methods have gained popularity recently and has been used to study problems such as modeling phospholipidic vesicles and cardiomyocyte membrane (Maitre et al., 2009), modeling large number of red blood cells (RBCs) and platelets in a capillary vessel (Sugiyama et al., 2017) and turbulent channel flow over hyperelastic walls (Rosti and Brandt, 2017). A recent work by Kamrin et al. (2012) introduced the “reference map technique” (RMT), a fully Eulerian approach for the simulation of solids and an extension to coupled fluid-solid problems (Valkov et al., 2015). In this work, visco-elastic solids were successfully simulated on a staggered grid coupled with a Newtonian fluid in a compressible flow setting using hyperelastic constitutive laws. The main differences of the RMT and VOF methods are (i) a reference-map vector field is transported to track the deformation of the solid in RMT approach as opposed to a tensor field (left Cauchy-Green deformation tensor) in the VOF method. (ii) RMT method has been extended to account for solid-solid contact conditions. We therefore adopt the RMT approach and extend this formulation for incompressible settings (Jain and Mani, 2017; Rycroft et al., 2018) and assess multiple improvements to the original RMT method (Valkov et al., 2015). Other approaches in the literature similar to RMT that is worth mentioning are (Dunne, 2006; Govindjee and Mihalic, 1996). Further, we point the readers towards excellent review articles by Hou et al. (2012); Takizawa et al. (2012); Takagi et al. (2012); Mittal and Iaccarino (2005) for additional details on the methods and many more applications of FSI problems which could not be included here for the sake of brevity.

In the present paper, we describe a conservative and non-dissipative reference-map technique (RMT) for the simulation of incompressible soft solids in fluids. We discuss the improvements made for this model in terms of the accuracy, cost, ease of implementation, and robustness of the method and also discuss some of the best modeling practices. Some of the important features of our approach compared to the state-of-the-art RMT (Valkov et al., 2015) are (a) discrete momentum conservation, (b) a least-squares extrapolation procedure that is accurate and cost-effective, (c) a modified advection equation for the reference map field that improves robustness of the method, (d) a non-dissipative central-difference scheme that eliminates any spurious dissipation of kinetic energy, and (d) projection method for incompressible flows. Rest of the paper is organized into sections as follows: Section 2 describes the basic formulation of the reference-map technique, governing equations that describe the motion of fluids and solids, and their respective constitutive laws. Section 3 describes the numerical method and introduces the conservative formulation, discretizations, projection method algorithm, a strategy to reconstruct level-set field, modifications to the reference-map advection equation, a new least-squares based extrapolation procedure and a closure model. Section 4 presents the verification of the solver against the results from a Lagrangian approach, presents the cost and accuracy improvements of the new extrapolation procedure, illustrates the importance of the use of a conservative formulation, and presents more complex test cases involving solid-solid and solid-wall contact situations. Finally, section 5 presents the summary along with the concluding remarks.

2 Eulerian formulation for solids and fluids

2.1 Reference map technique

Consider a solid in convective coordinate system, as shown in Figure 2. At time $t=0$, the solid is in its initial configuration (reference configuration), represented by $\Omega_{0}$, and at time $t>0$, the solid is in its deformed configuration, represented by $\Omega_{t}$ after being displaced and deformed by external forces. If, $\vec{X}$ represents a position vector in $\Omega_{0}$ that points to a material particle, then this material particle in $\Omega_{t}$ has the same $\vec{X}$ associated with it, since $\vec{X}$ represents the initial coordinates of the point in $\Omega_{0}$. Hence $\vec{X}$ acts as a tag for all the material particles in the solid. If $\vec{x}$ represents the corresponding position vector in $\Omega_{t}$, then we can define a vector map $\vec{\xi}:\mathbb{R}^{4}\rightarrow\mathbb{R}^{3}$ (a reference map) as

[TABLE]

such that $\vec{\xi}$ remains constant for a material particle in the solid (as long as the solid doesn’t deform plastically) but varies from particle to particle. Hence the material derivative of $\vec{\xi}$ field yields

[TABLE]

Expressing this in terms of the local derivatives, we obtain an advection equation for the $\vec{\xi}(\vec{x},t)$ field as

[TABLE]

This equation can be integrated in time given the initial condition $\vec{\xi}(\vec{x},t=0)=\vec{x}=\vec{X}$. Thus, $\vec{\xi}(\vec{x},t)$ acts as a tag for all the points in the solid, and the kinematic condition in Eq. (3) can be used to track every point in the solid, given its initial coordinates. Stress and strain in solid constitutive laws are typically expressed in terms of the material deformation gradient $\mathbb{F}$. Hence, relating $\mathbb{F}$ to $\vec{\xi}(\vec{x},t)$ as

[TABLE]

we can express the stress and strain tensors in terms of this new primitive variable $\vec{\xi}(\vec{x},t)$. Eqs. (1)-(4) in combination give rise to a novel approach to track all the material points in a solid and close the system of equations to model a solid on an Eulerian grid.

2.2 Governing equations for solids and fluids

In an Eulerian formulation, momentum balance equation for both fluids and solids can be written as

[TABLE]

where $\vec{u}$ is the global velocity field and $\underline{\underline{\sigma}}$ is the Cauchy stress. Mass balance equation for fluids (continuity equation) can be written as

[TABLE]

which in the incompressible limit simplifies to $\vec{\nabla}.\vec{u}=0$. Similarly, the mass balance for solids can be written as $\rho=\rho_{0}[det(\mathbb{F})]^{-1}$, and in the incompressible limit it simplifies to $det(\mathbb{F})=1$, implying that the density doesn’t change ($\rho=\rho_{0}$). Here, $\rho$ and $\rho_{o}$ are the density of the deformed and reference configurations, respectively. It can be shown that the conditions $\vec{\nabla}.\vec{u}=0$ and $det(\mathbb{F})=1$ are equivalent (see, Appendix A).

For solids, the Cauchy stress can be expressed as a function of strain given by

[TABLE]

where $\mathbb{E}=(1/2)(\mathbb{F}^{T}\mathbb{F}-\mathds{1})$ is the Green’s (or Lagrangian) finite strain tensor, $\mathbb{C}=\mathbb{F}^{T}\mathbb{F}$ is the right Cauchy-Green’s deformation tensor (or stretch tensor), $\psi(\mathbb{F})=\bar{\psi}(\mathbb{E})=\hat{\psi}(\mathbb{C})$ is the strain energy density (Helmhotz free-energy density) function and $\lambda=P$ is the Lagrangian multiplier and is equal to pressure in the incompressible limit (Holzapfel, 2000). We use the incompressible neo-Hookean constitutive model for solids, given by

[TABLE]

where $\mu^{s}={E}/{2(1+\nu)}$ is the shear modulus (Lame’s first parameter), $E$ is the Young’s modulus and $\nu$ is the Poisson’s ratio. Taking a partial derivative of this strain energy density function with respect to $\mathbb{C}$ yields

[TABLE]

Using this and the incompressibility condition for solids ($det(\mathbb{F})=1$), the $\underline{\underline{\sigma}}^{s}$ reduces to a simple form given by $\underline{\underline{\sigma}}^{s}=2\mu^{s}\mathbbm{b}-P\mathds{1}$, where $\mathbbm{b}=\mathbb{F}\mathbb{F}^{T}$ is the left Cauchy-Green’s deformation tensor (or stretch tensor). Further more, expressing $\mathbb{F}$ in terms of $\vec{\xi}$, Cauchy stress can be expressed in terms of this new primitive variable $\vec{\xi}$ as

[TABLE]

Nonlinearity in the stress-strain relationship is more evident when $\vec{\xi}$ is expressed in terms of its components. For an incompressible solid in two-dimensions, the Cauchy stress reduces to the form (Appendix B)

[TABLE]

where $\alpha=\vec{\xi}.\hat{i}$ and $\beta=\vec{\xi}.\hat{j}$ are the components of $\vec{\xi}$. For fluids, the Cauchy stress can be expressed as a function of the rate of strain. We use the Newtonian constitutive model given by

[TABLE]

where $P$ is the pressure and the matrix form of the above system of equations is for a fluid in two-dimensions. Extension to three-dimensions is not included here, but is straightforward. We solve a conservative variable-density formulation of the above system of equations and to close the system of equations for fluid-solid coupled simulations, we use a mixture model derived based on the one-fluid formulation (Kataoka, 1986) for two-phase flows (see, Section 3.6), or the so called “one-continuum formulation” (Sugiyama et al., 2011).

3 Discretization, numerical method and conservative implementation

3.1 Basic methodology

Consider a solid on an Eulerian grid, as shown in Figure 3. Here, $\Omega_{S}$ represents the region inside the solid, $\partial\Omega_{S}$ represents the boundary of the solid, $\pm\Omega_{E}\approx 4\Delta x$ represents a narrow band region (an extended solid region) around $\partial\Omega_{S}$, and $\pm\Omega_{T}\approx 2\Delta x$ represents another narrow band region (a transition zone) around $\partial\Omega_{S}$ (not shown in the Figure 3). Since, both solid and fluid regions are solved together in a coupled fashion, they share the same grid and a global velocity field. In the regions of solid $\Omega_{S}$, solid Cauchy stress $\underline{\underline{\sigma}}^{s}$ is computed using the solid constitutive law (Eq. 11), and outside this region, fluid Cauchy stress $\underline{\underline{\sigma}}^{f}$ is computed using the fluid constitutive law (Eq. 12). Once the stresses for solid and fluid regions are evaluated, a level-set field $\phi$ and a smoothed Heaviside function $H(x)$ is constructed using the reference map field $\vec{\xi}$ as illustrated in Figure 4. This Heaviside function is used to appropriately blend the solid and fluid stresses around the solid-fluid interface to compute the global Cauchy stress $\underline{\underline{\sigma}}$. Finally the velocity field is updated by solving the discretized version of momentum equation and by projecting the velocity field onto a divergence-free field. Form of the equations used, discretization techniques and algorithms used in this approach are explained in detail in the subsequent sections.

An important thing to note is that $\vec{\xi}$ is a variable that contains the information of the origin of the material. This field quickly becomes invalid in the region containing fluids due to the highly nonlinear deformation behavior of the fluids. Hence, $\vec{\xi}$ is defined only within the solid region, and to evaluate the solid stress in $\Omega_{T}$, the $\vec{\xi}$ field is appropriately extrapolated into the regions outside the solid (into $\Omega_{E}$).

Valkov et al. (2015), used the hyperbolic partial-differential equation (PDE) approach of Aslam (2004) to extrapolate the $\vec{\xi}$ field. This approach assumes that a level-set field is known in the region of extrapolation. By contrast, we use a least-squares-based extrapolation procedure (see, Section 3.5) that does not require a known level-set field in the region of extrapolation. However, a local level-set field $\phi$ has to be defined in the $\Omega_{E}$ region, which is used in defining the Heaviside function required for the mixture model (see, Section 3.6) and also in enforcing the solid-solid and solid-wall contact boundaries. Hereafter, we refer to the approach by Valkov et al. (2015) as the original RMT.

3.2 Conservative formulation and discretization

In the numerical solution of partial differential equations, divergence form of the equations is usually preferred over the primitive form (non-conservative form), since it results in discrete conservation of the quantities being solved. We solve the momentum equation in a conservative form as written in Eq. 5, where both the inertial term and the stress term are in divergence form (see, Section 4.3 for the illustration of importance of the use of divergence form for the stress term). We also use the conservative form of the equation for the advection of $\vec{\xi}$ field. Though $\vec{\xi}$ is not a physically conservative field, volume enclosed in the solid region $\Omega_{S}$ bounded by the fluid-solid interface $\partial\Omega_{S}$ that is extracted using the $\vec{\xi}$ field should be conserved (see, Section 4.6.1 for a description on the volumetric error of the solid). Notice that the Eq. 3 can be rewritten in a conservative form as

[TABLE]

in the incompressible limit, using the divergence-free condition ($\vec{\nabla}\cdot\vec{u}=0$).

We use a finite-volume approach on a collocated uniform grid to discretize our system of equations. Hence, all our primary variables ($\vec{\xi}$,$\vec{u}$,$p$,$\rho$) are stored on the cell center. We modify the approximate projection method of Almgren et al. (2000) to incorporate the coupled solution of solid and fluid regions. The steps involved in our projection method are shown in detail in Algorithm 1.

We split the momentum equation into an advection and diffusion part, since it allows us to use different time-stepping schemes. A second-order central differencing scheme is used to compute convective fluxes in the advection part of the momentum equation and the advection equation for $\vec{\xi}$, and they are solved using an RK4 time integration scheme. The use of central-difference scheme for the advection of both $\vec{\xi}$ and $\vec{u}$ fields not only yields a conservative and a non-dissipative approach but also results in solving momentum equation and reference map advection equations consistently, which is crucial for the simulation of high solid-to-fluid density-ratio flows (see Section 3.10 in Tryggvason et al. (2011)). A forward-Euler time integration scheme is used to solve the diffusion part of the momentum equation (Eq. 15). We use the second-order central-difference approximation to evaluate the gradient tensor ($\vec{\nabla}\vec{\xi}$) in Eq. (10), unlike the one-sided differences used in the original RMT. For example, in a Cartesian two-dimensional case, $\partial\alpha/\partial x$ in Eq. (11) is approximated as

[TABLE]

where $\alpha=\vec{\xi}.\hat{i}$. The divergence of Cauchy stress in Eq. (15) is also computed using the second-order central-difference scheme. The exact form of discretization of the stress terms is crucial in obtaining a consistent and conservative formulation, hence the discretization used in the current work is presented in detail in Appendix C.

Finally, the use of collocated grid arrangement results in checkerboard pressure fields. To eliminate this, we use a “Rhie-Chow like interpolation” for the intermediate velocity fields ($\vec{u}^{**}$) after the update of advection and diffusion

[TABLE]

where subscript $P$ represents cell-centered values, $f$ represents face-centered values, $\left\langle\right\rangle_{P\rightarrow f}$ is an interpolation from cell center to the cell face and $F$ is the body force computed using the balanced-force approach of Francois et al. (2006)). Interpolation from cell center to the cell face ($\left\langle\right\rangle_{P\rightarrow f}$) for a uniform Cartesian grid can be written as

[TABLE]

[TABLE]

where $u$ and $v$ are the $x$ and $y$ components of the velocity field. This Rhie-Chow like interpolation does not affect the discrete conservation of momentum. However, it does add a small amount of conservation error in the transport of kinetic energy which is of the order $O(\Delta t\Delta x^{2})$. This has been previously shown to be of dissipative in nature, hence it does not affect the stability of the method (see, Section 6.1 of Morinishi et al., 1998).

3.3 Level-set reconstruction

As explained in the section 3.1, a level-set field $\phi$ is required to be defined at every time step in the $\Omega_{E}$ region. One way to define $\phi(\vec{x},t)$ is to advect $\phi$ using the standard level-set advection equations given $\phi(\vec{X},t=0)$. This approach could lead to a mismatch between the $\phi(\vec{x},t)=0$ and the boundary of the solid defined by $\vec{\xi}(\vec{x},t)$ field, which in turn could result in a wrinkled solid-fluid interface, affecting the overall quality of $\vec{\xi}$ field in the extrapolated region (see Figures 9,10 in Valkov et al. (2015)). In the original RMT approach, this issue was resolved by performing additional smoothing routines to eliminate the striations in the extrapolated regions, which could potentially lead to additional mass conservation issues. To avoid this problem, we propose a simpler, exact, conservative and also cost-effective way to define the level-set field $\phi(\vec{x},t)$ at any time $t$. $\phi(\vec{x},t)$ can be reconstructed from the given $\phi(\vec{X},0)$ field at $t=0$, utilizing the known $\vec{\xi}(\vec{x},t)$ field at time $t$ using a simple condition given by

[TABLE]

Since an analytical expression can be defined for $\phi(\vec{X},t=0)$ for simple-shaped solids, the above equation yields an exact field for $\phi(\vec{x},t)$ for a given $\vec{\xi}(\vec{x},t)$, thus maintaining a perfect match between the $\phi(\vec{x},t)=0$ surface and the boundary of the solid defined by $\vec{\xi}(\vec{x},t)$, which is crucial in developing a robust solver. If an analytical expression for $\phi(\vec{X},t=0)$ is not available, then a bilinear interpolation (in two-dimensions) can be used to calculate $\phi(\vec{\xi}(\vec{x},t),t=0)$.

3.4 Modified reference map advection

The reference map field $\vec{\xi}$ is advected using Eq. (13). As explained in Section 3.1, $\vec{\xi}$ is defined and advected only within the $\Omega_{S}$. This can be conveniently achieved by modifying Eq. (13) into

[TABLE]

where $H(\vec{x})$ is a Heaviside function defined as

[TABLE]

This modification to the advection equation of $\vec{\xi}$ has multiple advantages; the very obvious one is that this approach effectively eliminates the high-frequency content in the $\vec{\xi}$ field, resulting in an ability to use simple schemes such as central-differences to compute the fluxes, without losing the accuracy of the solution to dispersion errors. To realize the second advantage, which is more subtle, consider the $\vec{\xi}$ field of a one-dimensional solid, as shown in Figure 5.

At time $t=0$, $\vec{\xi}$ is a simple straight line ($\xi=x$) and is given as an input to the solver, as shown on the left. Let $u$ denote the velocity field; then after time $t$, an ideal solid would have advected to a new location, shown on the right, maintaining the shape (solid line). If the modified equation shown in Eq. (26) is not used to advect, then a Total Variation Diminishing (TVD) type scheme can be used to compute the fluxes, which artificially add diffusivity to the equation to stabilize the solver. This results in a non-monotonic $\vec{\xi}$ profile shown on the right with dotted lines. If this profile is obtained as a result of advection, then the reconstruction step of the level-set field using Eq. (25) breaks down. To understand this, consider two solid circles 1 and 2 in Figure 5. Circle 1 is inside the solid in $\Omega_{S}$ at $t=0$, but circle 2 is outside the solid in $\Omega_{E}$. After the advection step, circle 1 still represents a value of $\xi$ inside the solid region, whereas circle 2 now represents a value of $\xi$ inside solid region $\Omega_{S}$ due to numerical diffusion. To clip the values of $\xi$ outside the solid region before extrapolating $\xi$, the boundaries of the solid needs to be identified. This can be done using the level-set field constructed using Eq. (25) (which takes $\xi$ as the input). This procedure (without the modification to the advection equation Eq. 26) typically creates two boundaries for the solid, resulting in the failure of the method. Therefore using the modified advection equation for $\vec{\xi}$ effectively eliminates this issue by clipping the values of $\vec{\xi}$ outside the solid right in the advection step. As a result, $\vec{\xi}$ can be extrapolated without any need for explicit clipping.

3.5 Least-squares-based extrapolation

The original RMT used a hyperbolic PDE approach to extrapolate the $\vec{\xi}$ field outside the solid regions into $\Omega_{E}$. We propose a simpler, more cost-effective approach to extrapolate the $\vec{\xi}$ field based on the assumption that the $\vec{\xi}$ field is locally linear, even in the deformed state of the solid since the $\vec{\xi}$ field represents a mapping from the current coordinates $\vec{x}$ to the original coordinates $\vec{X}$ of the solid, which is linear in $x$, $y$ and $z$. This assumption relies on the fact that the solids do not continuously deform under an applied stress. Unlike solids, liquids do not resist stress and continuously deform which would violate this locally linear assumption of $\vec{\xi}$ field. Therefore the $\vec{\xi}$ field is only defined within the solid region. This assumption of locally linear $\vec{\xi}$ field was also used in the hyperbolic PDE based extrapolation of the original RMT approach (Valkov et al., 2015).

Consider the solid represented by a square (in two-dimensions) in Figure 6. Consider the dashed circle of radius 4$r$ as the stencil, where $r=\sqrt[]{\Delta x^{2}+\Delta y^{2}}$. Hence a plane of the form $\xi=ax+by+c$ can be fit for the known cell values, where $x$, $y$ and $\xi$ represent the coordinate location and the reference map value of the cells and $a$, $b$ and $c$ are the coefficients to be determined, thus forming an over-determined system that can be solved using the least-squares approach. The stencil’s radius was chosen to make the system over-determined for all the possible configurations. Once the coefficients are calculated, the value of $\xi$ at the solid circle can be computed. The procedure begins by repeatedly solving least-squares systems for all the cells adjacent to the cells for which the value of $\xi$ is already known. This is considered as the first pass. The values computed in the first pass are considered as good as the values inside the solid for the second pass. This procedure is repeated until the required width of the extrapolated region is obtained. This cell traversal procedure is summarized in Algorithm 2.

3.6 Closure model

The level-set field $\phi$ reconstructed using Eq. (25) should be reinitialized to restore its signed-distance property. We solve the Eikonal equation by adopting the fast marching method (FMM) of Chopp (2001) to reinitialize the $\phi$ field. The coupled fluid-solid system of equations is closed by defining the mixture model inspired by the “one-fluid formulation” as

[TABLE]

[TABLE]

where $\hat{\phi}$ is the reinitialized level-set field and $\hat{H}(x)$ represents a smoothed Heaviside function defined as

[TABLE]

where $w_{T}$ represents the width of the transition region $\Omega_{T}$. For $n$ number of solids, this model can be extended accordingly

[TABLE]

When two solids collide in a fluid, a body force needs to be added to the momentum equation to keep them separated and to avoid the inter-penetration of solids. We use a similar procedure as described in Valkov et al. (2015) to calculate the body force $\vec{f}_{i,j}$ for solid-solid contact and solid-wall contact conditions. A level-set field $\phi_{12}$ is defined as

[TABLE]

where $\phi_{1}$ and $\phi_{2}$ are the level-set fields associated with two colliding solids, hence $\phi_{12}=0$ represents a mid-surface between the two solids. The body force $\vec{f}_{i,j}$ can then be defined as

[TABLE]

[TABLE]

where $\hat{n}_{12i,j}$ is the unit vector normal to the level-sets of $\phi_{12}$ and pointing away from the mid-surface, $k_{rep}$ is a prefactor and $\delta_{s}(x)$ is a compactly supported influence function given by

[TABLE]

Note that, we need to define a separate $\vec{\xi}$ field and transport it for each object that undergoes collision in the simulation, which is required to evaluate Eq. (30). However, if there are many objects in the domain and if we know that they do not collide with each other beforehand (when they are sufficiently far away), we can use the same $\vec{\xi}$ field for them to reduce the cost and memory requirements.

Finally, the pressure Poisson equation (Eq. (17)) which results in a linear system of equations is solved using a conjugate-gradient (CG) approach. Major differences and improvements to the original RMT method by Valkov et al. (2015) are listed in Table 1. We thus extended the original reference map technique (RMT) to solve for incompressible fluid-structure interaction problems on an Eulerian collocated grid. Modifications proposed in the extrapolation procedure of the reference map, reconstruction of the level-set field and consistent numerical discretization results in improved robustness, cost effectiveness and conservation properties of the approach.

4 Results and discussion

In this section, we first present some basic validation test cases to assess the accuracy and cost of our fluid and coupled fluid-solid solver. This is then followed by more complex cases involving solid-solid and fluid-solid contact conditions. Since the fluid-solid coupled problem involves multiple time scales, and an explicit time integration is adopted for solving the system of equations, care must be taken to satisfy all the time step constraints involved in the problem. Time step restriction due to CFL criterion from the advection can be written (for forward Euler in one dimension) as $\Delta t\leq C\Delta x/u$, where $C$ represents the Courant number. Time step restriction from the diffusion equation for fluids yields $\Delta t\leq 0.5\rho^{f}(\Delta x)^{2}/\mu^{f}$. Similarly, shear waves in the solids need to be resolved, and the speed of this shear wave is given by, $u=\sqrt[]{\mu^{s}/\rho^{s}}$. Hence a time constraint based on this shear wave speed can be defined as $\Delta t\leq P\Delta x\sqrt[]{\rho^{s}/\mu^{s}}$, where $P$ represents an appropriate pre-factor that depends on the numerical method. If the ratio of $\mu^{s}/\rho^{s}$ is high, travelling shear waves in the solid typically imposes the most restrictive time constraint of all. Hence, in the stiff solid limit such as in the metals, imposed time step constraints are so strict that the simulation time close to solid length scales is virtually impossible with the explicit time stepping approach. Therefore, this formulation is best suited for the simulation of soft solids in fluids.

4.1 Validation of the fluid solver

The incompressible Navier-Stokes solver on a collocated grid was validated for the lid-driven cavity case against the benchmark results from Ghia et al. (1982). A $100\times 100$ grid was used for the simulation, and the results are reported for $\mathrm{Re}=1000$. Figure 7 shows a good match of the $u$ and $v$ velocities along the vertical and horizontal lines through the center of the domain with the results from Ghia et al. (1982).

Since the equations are solved on a collocated grid, to eliminate the checkerboard fields a Rhie-Chow-like interpolation was performed, as described in Section 3.2. Figure 8 presents pseudocolor plots of the velocity and pressure fields from the lid-driven cavity case, illustrating the smoothness of the solution fields obtained.

4.2 Cost and accuracy of the extrapolation procedure

We compared the accuracy of our least-squares extrapolation procedure with that of the hyperbolic partial differential equation (PDE) approach used in RMT (using a RK2-minmod scheme to solve the hyperbolic PDEs). Figure 9 shows the results of the Zalesak disk test case, wherein a slotted disk (a $\vec{\xi}$ field) that is placed off-center is advected with a given background rotational velocity field and compared against the initial conditions after one full rotation. Three solid lines in (a) represent the initial and final $\phi$ fields of $\partial\Omega_{S}$ (fluid-solid interface), $+\partial\Omega_{E}$ and $-\partial\Omega_{E}$ (boundaries of the extended solid region). Solid lines in (b) and (c) represents $\partial\Omega_{S}$ and the shaded region represents $\alpha=\vec{\xi}\cdot\hat{i}$ and $\beta=\vec{\xi}\cdot\hat{j}$ fields.

The initial and final $\phi$ fields in Figure 9 (a) are exactly on top of each other, showing that the extrapolation procedure by itself is very accurate. One should be careful in interpreting this result, and should not relate this with the rotation of Zalesak disk usually presented in the literature that is obtained as a result of direct advection of $\phi$ field . Here $\phi$ field is reconstructed using the condition in the Eq. 25 and the high accuracy of this $\phi$ field could only be achieved due to the advection of $\vec{\xi}$ that was linear and smooth using a second-order central scheme. This clearly shows the advantage of using the compatibility condition in the Eq. 25 as opposed to the advection of the $\phi$ field.

Further, since the errors in the extrapolation procedure manifests as the error in the advection of the $\vec{\xi}$ field, we computed the $L_{2}$ norm error $E_{\xi}=||\vec{\xi_{i}}-\vec{\xi_{f}}||_{2}$ for the advection, where $\vec{\xi_{i}}$ and $\vec{\xi_{f}}$ are the initial and final fields obtained after one full rotation, and report them in Table 2. It is evident that our least-squares procedure is considerably more accurate when compared to the PDE approach.

The above test case was performed on a $100\times 100$ grid. Moreover, we also compared the cost of the extrapolation procedure using both the approaches and found that on an average the least-squares procedure required $\approx 100\mathrm{ms}$ per extrapolation, whereas the PDE approach required $\approx 1550\mathrm{ms}$ per extrapolation on this grid (close to the time taken by a Poisson solver), for an extrapolation band region of $5\Delta x$. This also proves that our least-squares procedure is extremely cost-effective when compared to the PDE approach.

Additionally, to demonstrate that the use of compatibility condition in the Eq. 25 to reconstruct $\phi$ field is not limited to simple shapes that have analytical expression, we considered an asymmetric star-looking object that has sharp regions and repeated the exercise above. Figure 10 shows the results of the one full rotation of the star-looking object advected with a given background rotational velocity field. Three solid lines represent the initial and final $\phi$ fields of $\partial\Omega_{S}$ (fluid-solid interface), $+\partial\Omega_{E}$ and $-\partial\Omega_{E}$ (boundaries of the extended solid region). The initial and final $\phi$ fields in Figure 10 (a) are again exactly on top of each other, showing the high accuracy of the method even for objects with sharp regions. Further to quantify the error, we computed the $L_{2}$ norm error $E_{\xi}=||\vec{\xi_{i}}-\vec{\xi_{f}}||_{2}$, where $\vec{\xi_{i}}$ and $\vec{\xi_{f}}$ are the initial and final fields obtained after one full rotation. The error values are $E_{\vec{\xi}.\hat{i}}=6.39\times 10^{-9}$ and $E_{\vec{\xi}.\hat{i}}=6.77\times 10^{-9}$, which are of the same order as the ones reported in Table 2 for a grid of size $100\times 100$.

Finally to study the effect of grid size on the sharp corners, we repeated the same test case for various grids of sizes $50\times 50,100\times 100,200\times 200$ and $400\times 400$. The final shape of the object for various grid sizes is shown in Figure 11, along with close-up views around a sharp corner and a smooth corner showing the grid convergence. We also computed a normalized $L_{2}$ norm error $E_{\xi norm}=||\vec{\xi_{i}}-\vec{\xi_{f}}||_{2}/(N_{x}\times N_{y})$, where $N_{x}$ and $N_{y}$ are the number of grid points along $x$ and $y$ directions, and report them in Table 3. Normalization is done in such a way that the error quantity $E_{\xi norm}$ being compared is grid-size independent and that it represents the error incurred per grid cell in the domain. Clearly, the error per grid cell is very close to machine accuracy for all grid sizes and are roughly independent of the grid size.

4.3 Conservative vs Non-conservative implementation

Here we would like to highlight that a careful implementation of the blending of fluid and solid Cauchy stresses is crucial in obtaining a discretely conservative momentum formulation. For example, one approach is to compute fluid and solid Cauchy stresses ($\underline{\underline{\sigma}}^{s},\underline{\underline{\sigma}}^{f}$), combine them to obtain a global Cauchy stress ($\underline{\underline{\sigma}}$), and then calculate the divergence of this stress to obtain the force per unit volume ($\vec{f}$) due to stresses as

[TABLE]

[TABLE]

The second approach is to compute the divergence of the solid and fluid Cauchy stresses ($\vec{\nabla}\cdot\underline{\underline{\sigma}}^{s},\vec{\nabla}\cdot\underline{\underline{\sigma}}^{f}$) and combine them to obtain the force per unit volume as

[TABLE]

The first approach is the one that leads to a conservative formulation, due to the presence of divergence outside the blending operation. This divergence operator, when summed up over adjacent control volumes, leads to an exact cancellation of the terms (analogous to a telescoping series). Hence, we use the conservative formulation in our solver.

A simulation of a solid placed in a Taylor-Green vortex was performed to qualitatively study the differences between these two formulations. Consider Figure 12, which shows the initial state of a solid placed in a Taylor-Green vortex field. Initial flow field should stretch the solid to a certain extent, beyond which the internal stresses developed in the solid should retract it back resulting in an oscillating motion of the solid that stretches and retracts back and forth until all the energy is lost in the viscous dissipation of the fluid. Figure 13 shows the result of the simulation performed using both the non-conservative formulation and conservative formulation described above. Clearly, the results are completely unphysical for the non-conservative formulation wherein the solid extends indefinitely with no signs of retraction. By contrast, the conservative formulation for the exact same problem resulted in a more physically meaningful calculation. This simple demonstration illustrates the importance of a conservative numerical implementation, very much similar to the one in compressible flows to achieve correct shock speeds (Laney, 1998) and in high-density ratio two-phase flows (see Figure 7 in Raessi and Pitsch (2012)).

4.4 Convergence study

Above demonstrated test case of a solid placed in a Taylor-Green vortex field was repeated for the values used in (Zhao et al., 2008; Robinson-Mosher et al., 2011) to validate our solver against the results from a mixed Eulerian-Lagrangian based approach. A solid of radius $r=0.2$ is placed in an initially imposed Taylor-Green vortex field given by the streamfunction $\psi=\psi_{0}sin(k_{x}x)sin(k_{y}y)$ where $\psi_{0}=5\times 10^{-2}$ and $k_{x}=k_{y}=2\pi$. Domain size used is $1\times 1$ and is discretized into a $128\times 128$ grid. Other parameters used in the simulation are fluid viscosity $\mu^{f}=10^{-3}$, shear modulus $\mu^{s}=0.5$, solid density $\rho^{s}=1$ and fluid density $\rho^{f}=1$. For the sake of consistency with the results of Robinson-Mosher et al. (2011); Zhao et al. (2008), a small amount of viscosity equal to the fluid viscosity of $\mu^{f}=10^{-3}$ is added in the solid regions. But in general our solver is stable without any viscous damping in the solid regions (see section 4.6 for simulations without any viscosity in the solid regions). Time evolution of kinetic energy ($ke$) and strain energy ($se$) is plotted in Figure 14 for various grid size and also compared against previous studies, where

[TABLE]

and

[TABLE]

Clearly the results are independent of the grid for sizes $128\times 128$ and above. Further, the viscous dissipation ($\varepsilon$) in the fluid and solid regions combined was computed using the expression

[TABLE]

and the conservation of total energy $E$ was assessed at the final time of $t=1$ and was found to decrease less than $1\%$ of the initial time value, where $E$ is given by

[TABLE]

Frequency of oscillation of the solid matches well with the results of Robinson-Mosher et al. (2011); Zhao et al. (2008). The time evolution of the kinetic energy in the present work matches well with that of Robinson-Mosher et al. (2011) during the first period of oscillation, but eventually the kinetic energy decays faster in the simulations by Robinson-Mosher et al. (2011); Zhao et al. (2008) compared to the current results. Similarly, the strain energy is under-predicted and decays faster in the previous works compared to the current results. This highlights the non-dissipative nature of central-difference scheme used in the current work. Further, using the same test case we also assessed the order of convergence of all the primitive variables ($\vec{u},p,\vec{\xi}$), the kinetic energy ($ke$) and the strain energy ($se$) of the solid used in our solver against a refined case on a $1024\times 1024$ grid. Figure 15 shows that the order of convergence is roughly $O(\Delta x^{2})$ for all the variables. Errors are defined as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the subscript $ref$ refers to the most refined case on a $1024\times 1024$ grid.

4.5 Solid in a driven cavity

To further validate the fluid-solid coupling of the solver, we simulated a deformable solid in a lid-driven cavity. This test case was previously simulated using mixed-Lagrangian-Eulerian based approach by Zhao et al. (2008) and using a VOF based Eulerian approach by Sugiyama et al. (2011). Figure 16(a) shows the initial configuration of the solid in the domain. Domain used for this simulation is $[0,1]\times[0,1]$ and is discretized into a $128\times 128$ grid. The solid is initially circular in shape with a radius of $r=0.2$ and is placed at ($0.6,0.5$) location. Other parameters used in the simulation are fluid viscosity $\mu^{f}=10^{-2}$, solid viscosity $10^{-2}$, shear modulus $\mu^{s}=0.05$, solid density $\rho^{s}=1$ and fluid density $\rho^{f}=1$. Time evolution of the interface of two solids are shown in Figure 16, where the black solid line represents the current method and the dashed red represents results by Sugiyama et al. (2011), which shows a pretty good match. Further, we also plot the centroid of the solid in space in Figure 17 against the results by Sugiyama et al. (2011). This shows that the centroid obtained using the present conservative Reference-Map-Technique on a grid of $128\times 128$ is very close to the centroid obtained using a VOF based Eulerian method of Sugiyama et al. (2011) on a grid of $1024\times 1024$.

4.6 Simulations of solids in a fluid

In this section we present the simulations of more complex configurations of incompressible solid(s) in a fluid domain such as solid-solid contact, solid-wall contact situations. First, a case of solid-solid contact is considered. Figure 18 shows a configuration of two solids placed in an initially imposed Taylor-Green vortex field given by the streamfunction $\psi=\psi_{0}sin(k_{x}x)sin(k_{y}y)$ where $\psi_{0}=1$ and $k_{x}=k_{y}=1$. Domain used for this simulation is $[-\pi,\pi]\times[-\pi,\pi]$ and is discretized into a $100\times 100$ grid. Two solids are initially circular in shape with radii $r_{1}=r_{2}=\pi/3$ and are placed at ($\pi,1.4\pi$) and ($\pi,0.6\pi$) locations respectively. Other parameters used in the simulation are fluid viscosity $\mu^{f}=1$, shear modulii $\mu_{1}^{s}=\mu_{2}^{s}=100$, solid densities $\rho_{1}^{s}=\rho_{2}^{s}=100$ and fluid density $\rho^{f}=100$. Time evolution of the interface of two solids are shown in Figure 19. Solids collide and subsequently rebounce due to the internal stresses developed in them as a result of deformation. Centroid of both the solids are also plotted as a function of time.

The deformed configuration of the solids along with the normal stresses $\underline{\underline{\sigma}}_{n}=(\underline{\underline{\sigma}}_{11}+\underline{\underline{\sigma}}_{22})/2$ are shown in Figure 18 for the time $t=0.024$. Since the solids are in the rebouncing stage, formation of the four symmetric counter-rotating vortices can be clearly seen around the solids. A zoomed-in view of the solid is also shown in this Figure to illustrate the smoothness of the interface obtained in our approach even at such coarse resolution of $100\times 100$ grid points (due to the exact match between the level-set field $\phi$ and $\vec{\xi}$ fields at all times; see Section 3.3). We also do not see any striations in the extrapolated $\vec{\xi}$ fields that was observed in the original RMT (see Figure 10 in Valkov et al. (2015)), thus eliminating the requirement of the artificial smoothing routines that were used to remove the striations in the extrapolated region. This test case shows the robustness of our solver in handling the solid-solid contact situations. Further, the time evolution of the centroid of the colliding solids are plotted in Figure 20 for the grid sizes $100\times 100$, $200\times 200$ and $400\times 400$, which shows that the results (including the collision model) converge with the increase in grid size.

4.6.1 Deviation of $det(\mathbb{F})$ from $1$ (volumetric error)

In compressible flows, the conservative form of the momentum equation results in inconsistency between density advected using Eq. 6 and density computed using $\rho_{0}[det(\mathbb{F})]^{-1}$. To alleviate this, Kamrin et al. (2012) proposed an alternative approach for the computation of density in a one-dimensional setting, where the density is always defined in terms of the motion as opposed to solving the continuity equation. However, in the incompressible limit, density is constant within the solid. Therefore the Eq. 6 reduces to $\vec{\nabla}.\vec{u}=0$, which is satisfied discretely using the projection method. Hence the inconsistency is only in maintaining $\rho=\rho_{0}$ within the solid region, i.e., $det(\mathbb{F})$ equal to $1$. This condition is satisfied in the continous limit (see, Appendix A), and is generally not satisfied discretely. However, a good numerical implementation holds the value of $det(\mathbb{F})$ close to $1$. From Eq. (53) we can write

[TABLE]

Hence $det(\mathbb{F})-1=(dv-dV)/dV$ represents the local volumetric error in an incompressible solid due to the numerical discretization. Here, we present the volumetric error in the solid obtained using the present approach for the test case of two solids placed in an initially Taylor-Green vortex, described in Section 4.6. Figure 21 shows the error for the top solid at two different time instances (t=2.4, 5.6). Evidently, the volumetric error in the solid at $t=5.4$ is lower compared to that at $t=2.4$. Hence, the error ($det(\mathbb{F})-1$) does not seem to be accumulating with time, instead it is roughly proportional to the deformation of the solid. Moreover, the error is localized within the transition zone ($\Omega_{T}$) of the solid and the max value is around $0.14(14\%)$ and is independent of the grid size chosen. This error occurs due to the presence of mixture region where the stress is computed as a weighted average of the fluid and solid stresses and is typical of any Eulerian approach that uses a diffuse-interface approach. However, since this error is localized to the transition zone, the total error in the mean sense is negligible.

Additionally, to quantify the local volumetric error incurred throughout the solid the normalized $L_{2}$ norm of $det(\mathbb{F})$ from $1$ can be computed. This quantity is defined as $||det(\mathbb{F})_{i}-1||_{2}/n$, where $i$ is the cell index and $n$ is the number of cells inside the solid and is plotted as a function of time in Figure 22 (b) for three different grid sizes. Another similar measure that could be used to evaluate the deviation of $det(\mathbb{F})$ from $1$ is the net volumetric error of the solid, which represents the total volume loss or gain of the solid during the simulation. This quantity is defined as the normalized discrete summation of the $det(\mathbb{F})-1$ quantity, i.e., $\sum_{i=0}^{n}(det(\mathbb{F})_{i}-1)/n$ and is plotted as a function of time in Figure 22 (a). On a uniform grid this quantity can be expressed as

[TABLE]

where $V_{fin}$ and $V_{init}$ are the final and initial volumes of the solid. Figure 22 (a) shows the net volumetric error and Figure 22 (b) shows the normalized $L_{2}$ norm error as a function of simulation time for the top solid in the test case of two solids placed in an initially Taylor-Green vortex for various grid sizes. As the solid deforms, $det(\mathbb{F})$ deviates from $1$ and the volumetric error reaches a value of roughly $1\%$ and the $L_{2}$ norm error reaches a value of roughly $0.1\%$ for the $100\times 100$ grid case, however when the solid retracts back, $det(\mathbb{F})$ gets closer to $1$ and the volumetric error reduces down to $0.25\%$ and the $L_{2}$ norm error reduces to a value of roughly $0.02\%$. Therefore, there is no increase in error $det(\mathbb{F})-1$ or accumulation with time in the present approach. Furthermore, both volumetric and $L_{2}$ norm errors reduce significantly with increase in the grid size. Hence, the inconsistency does not pose a critical problem in the present method.

4.6.2 Simulations of solid-wall contact

Next, a sequence of three test cases named (a) collision, (b) bounce ($\mu^{s}=100$) and (c) bounce ($\mu^{s}=1000$) that involve solid-wall contact situations are considered. These classic test cases involving the collision of elastic solids with a rigid wall can be very useful and are of practical relevance in many engineering fields of research. In all the three cases, a domain of $[0,2\pi]\times[0,2\pi]$ is used and a circular shaped solid of radius $\pi/3$ is placed at ($\pi,\pi$) in an initially quiescent surrounding fluid. Case (a) is simulated in a microgravity condition ($g=0$) and the values of other simulation parameters used in this case are $\rho^{s}=100,\rho^{f}=100,\mu^{s}=100,\mu^{f}=1$. Solid is given an initial velocity of $\vec{u}=-1\hat{j}$ and since this initial condition is fictious and doesn’t satisfy incompressibility condition, the solver adjusts the velocity to achieve incompressibility in the first time step. Hence the effective velocity of the solid after one time step was $\vec{u}=0.48\hat{j}$. Time evolution of the interface of the solid is plotted as a function of time as shown in Figure 23. Solid encounters the rigid wall and bounces back and goes to a state of rest after losing all its kinetic energy to the surrounding fluid.

Case (b) is simulated in a gravity condition with $g=0.0981$ and the values of other simulation parameters used in this case are $\rho^{s}=1000,\rho^{f}=100$ hence a density ratio of $\rho^{s}/\rho^{f}=10$, $\mu^{s}=1000,\mu^{f}=10$. Solid is driven by the gravity and is initialized with a zero velocity. Time evolution of the interface of the solid is plotted as a function of time as shown in Figure 24. Solid encounters the rigid wall and bounces back and forth until it goes to a state of rest after losing all its kinetic and potential energy to the surrounding fluid. Case (c) is similar to case (b), but with parameters $\mu^{s}=100,\mu^{f}=1$. Solid is initialized with zero velocity and the time evolution of the interface of the solid is plotted as a function of time in Figure 25. Similar to case (b), here the solid bounces back and forth until it goes to a state of rest, but loses most of its energy to the fluid at its first encounter with the rigid wall due to a large deformation. The energy transferred to the fluid is eventually dissipated due to the action of viscosity. Though the strain energy stored in the solid in the event of a deformation is fully reversible/recoverable (non-viscous solid), energy spent in moving the surrounding fluid is large in the case of large deformations and hence the solid in the case (c), where the shear modulus is $\mu^{s}=100$, goes to rest much quicker when compared to the case (b) where the shear modulus is $\mu^{s}=1000$, with other parameters such as density ratio and gravity being identical. The centroid of the solid plotted as a function of the time for all three cases (a), (b) and (c) are shown in Figure 26.

5 Summary and Conclusions

We have presented an Eulerian formulation for the simulation of incompressible soft solids in a fluid. Methods that handle solids in a Lagrangian fashion are known to be too expensive for highly deforming solids due to large grid deformations and severe time step restrictions. On the other hand an Eulerian approach appears to be a more natural choice for such situations. Hence we have adopted the recently proposed “reference map technique” (RMT) by Valkov et al. (2015) to simulate solids and fluid-solid problems on an Eulerian grid. We extended this formulation for incompressible settings with the use of an approximate Projection method by (Almgren et al., 2000) to achieve divergence-free velocity condition.

Our formulation discretely conserves momentum and is very cost-effective. Furthermore, we introduced (a) a least-squares extrapolation procedure that is more accurate and cost-effective, (b) a modified advection equation for the reference map field that improves the robustness of the method, (c) a simple, cost-effective way to reconstruct the level-set field that removes any inconsistencies between the reference map field and the level-set field at all times and thereby eliminating the need to have more subroutines to fix the issue of striations of the interface (d) use of simple central-difference schemes to compute the fluxes that improves the stability of the numerical method and to eliminate any spurious dissipation of the kinetic energy.

We evaluated our solver on a variety of test cases involving solid-wall and solid-solid contact situations and showed that it is stable for all the cases. Furthermore, the test cases that we formulated can serve as a reference for future developers to compare and evaluate their models. Overall, this novel approach opens up a new pathway for the high fidelity numerical simulations of complex, large scale, coupled fluid-solid problems involving large deformations at lower costs compared to the Lagrangian or ALE methods.

Acknowledgments

This investigation was funded by the Office of Naval Research, Grant #119675. The first author was also supported by a Franklin P. and Caroline M. Johnson Fellowship. Authors would also like to thank the reviewer for her/his comments, that helped in improving this work.

Appendix A: Equivalence of $\vec{\nabla}\cdot\vec{u}=0$ and $det(\mathbb{F})=1$ relations

We can relate the divergence of velocity to the normal strain rate as

[TABLE]

where $\epsilon$ represents the strain rate. Therefore, the normal strain rate / total dilatation ($\epsilon_{ii}=0$) is zero because of the divergence free condition ($\vec{\nabla}\cdot\vec{u}=0$). Now, relating the total dilation to the rate of change of elementary volume with respect to a unit volume as

[TABLE]

we can therefore show that the rate of change of elementary volume ($\frac{D(\mathcal{V})}{Dt}$) is zero. This implies that the elementary volume $\mathcal{V}$ remains constant. Now rewriting this condition in a convective coordinate system, we obtain

[TABLE]

where $dv$ is the elementary volume after deformation and $dV$ is the elementary volume before deformation as shown in Figure 27. These two elementary volumes $dv$ and $dV$ can be re-expressed in terms of vectors $\vec{r}_{1}$, $\vec{r}_{2}$ and $\vec{r}_{3}$ that form the undeformed elementary volume $dV$ as

[TABLE]

[TABLE]

where $\vec{N}dA$ represents the area vector of the face of the elementary volume $dV$ formed by the vectors $\vec{r}_{1}$ and $\vec{r}_{2}$ and $\vec{n}da$ represents the area vector of the face of the elementary volume $dv$ formed by the vectors $\mathbb{F}\vec{r}_{1}$ and $\mathbb{F}\vec{r}_{2}$. Now making use of the Nanson’s formula ($\vec{n}da=det(\mathbb{F})\mathbb{F}^{-T}\vec{N}dA$) that relates these two area vectors, and making use of the relations in the Eqs. (51-52), we can show that

[TABLE]

Appendix B: Derivation of the incompressible solid Cauchy stress $\underline{\underline{\sigma}}^{s}$ in terms of components of the reference map $\vec{\xi}$

As described in Section 2.2, we can write the Cauchy stress in terms of the left Cauchy-Green’s deformation tensor (or stretch tensor, $\mathbbm{b}=\mathbb{F}\mathbb{F}^{T}$) for an incompressible neo-Hookean solid as

[TABLE]

Now, expressing $\mathbbm{b}$ in terms of $\vec{\xi}$ using the relation in Eq. 4, we obtain

[TABLE]

Further, rewriting $\vec{\xi}$ in terms of it’s components $\alpha=\vec{\xi}.\hat{i}$ and $\beta=\vec{\xi}.\hat{j}$ for a 2D system, we obtain

[TABLE]

which can be further simplified by evaluating the inverse of the matrix as

[TABLE]

We can further simplify the coefficient $C$ in the Eq. 56 using standard linear algebra identities and show that it is equal to $1$ with the use of the incompressibility condition for solids ($det(\mathbb{F})=1$) as

[TABLE]

Finally, substituting the expression for $\mathbbm{b}$ in Eq. 56 into the Eq. 54, we obtain the final expression for Cauchy stress $\underline{\underline{\sigma}}^{s}$ in terms of the components of the reference map field $\vec{\xi}$ as

[TABLE]

Appendix C: Discretization of the stress terms in the momentum equation

As shown in Appendix B, we can express the Cauchy stress in terms of the components of the reference map field $\vec{\xi}$ as

[TABLE]

where, $\alpha=\vec{\xi}.\hat{i}$ and $\beta=\vec{\xi}.\hat{j}$ for a 2D system. We use the second-order central difference scheme for the discretization of the gradients of $\alpha$ and $\beta$, which results in a conservative and non-dissipative formulation. If $i,j$ represents the cell index along $x$ and $y$ directions, we can write the discrete form of $\underline{\underline{\sigma}}^{s}$ as

[TABLE]

[TABLE]

[TABLE]

where $\sigma^{s}_{i,j,11}$, $\sigma^{s}_{i,j,12}$, $\sigma^{s}_{i,j,21}$ and $\sigma^{s}_{i,j,22}$ are the components of the tensor $\underline{\underline{\sigma}}^{s}_{i,j}$. Notice that we use a wider stencil that uses $i+1$ and $i-1$ points to obtain the gradient at $i$ as opposed to a compact stencil that uses $i+1/2$ and $i-1/2$. Once $\underline{\underline{\sigma}}^{s}_{i,j}$ is evaluated at $i,j$, we evaluate the fluid Cauchy stress $\underline{\underline{\sigma}}^{f}_{i,j}$ in a similar fashion using the same stencil. We then obtain the total Cauchy stress at $i,j$ as

[TABLE]

Further, to evaluate the divergence of $\underline{\underline{\sigma}}$ at $i,j$, we use the same stencil as

[TABLE]

which results in a conservative and consistent discretization of the stress terms that results in correct physical behavior of the solid as described in Section 4.3.