A moving control volume approach to computing hydrodynamic forces and torques on immersed bodies

TL;DR

This paper introduces a moving control volume method for accurately computing hydrodynamic forces and torques on complex bodies within an immersed boundary framework, improving robustness and applicability across flow regimes.

Contribution

It presents a novel moving control volume approach that simplifies force computation in immersed boundary methods, avoiding noisy derivatives and enabling use with adaptive meshes.

Findings

Method effectively suppresses force oscillations.

Excellent agreement with previous force/torque balance approaches.

Applicable across a range of flow regimes from Stokes to high Reynolds numbers.

Abstract

We present a moving control volume (CV) approach to computing hydrodynamic forces and torques on complex geometries. The method requires surface and volumetric integrals over a simple and regular Cartesian box that moves with an arbitrary velocity to enclose the body at all times. The moving box is aligned with Cartesian grid faces, which makes the integral evaluation straightforward in an immersed boundary (IB) framework. Discontinuous and noisy derivatives of velocity and pressure at the fluid-structure interface are avoided and far-field (smooth) velocity and pressure information is used. We re-visit the approach to compute hydrodynamic forces and torques through force/torque balance equation in a Lagrangian frame that some of us took in a prior work (Bhalla et al., J Comp Phys, 2013). We prove the equivalence of the two approaches for IB methods, thanks to the use of Peskin's delta functions. Both approaches are able to suppress spurious force oscillations and are in excellent agreement, as expected theoretically. Test cases ranging from Stokes to high Reynolds number regimes are considered. We discuss regridding issues for the moving CV method in an adaptive mesh refinement (AMR) context. The proposed moving CV method is not limited to a specific IB method and can also be used, for example, with embedded boundary methods.

Tables

Table 1. Table 1: Comparison of force measurements for Stokes flow around two concentric spheres.

Method	$\left|𝓕\cdot {𝐞}_x\right|$ (Inner)	$\left|𝓕\cdot {𝐞}_x\right|$ (Inner & Outer)
Analytical Eq. (36)	115.409	N/A
Lagrange multiplier Eq. (34)	114.297	17.350
Control volume Eq. (32)	114.298	17.350

Table 2. Table 2: Comparison of torque measurements 𝓜 ⋅ 𝐞 y ⋅ 𝓜 subscript 𝐞 𝑦 \bm{\mathrm{\mathcal{M}}}\cdot\bm{\mathrm{e}}_{y} for Stokes flow around single rotating sphere.

Grid size	Number of markers	$R_1^h/R_1^g$ [16]	Analytical Eq. (37)	LM Eq. (35)	CV Eq. (33)
${60}^3$	42	1.22	-269.669	-236.474	-236.474
${120}^3$	162	1.09	-192.323	-183.662	-183.662
${240}^3$	642	1.04	-167.052	-164.446	-164.446

Equations

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

$$\nabla \cdot \text{mbox}\mathrm{u}(\mathrm{x},t)$$

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

$$\mathrm{U}(\mathrm{s},t)$$

$$\frac{\partial \mathrm{X}}{\partial t}(\mathrm{s},t)$$

$$\rho \left(\frac{{\tilde{\mathrm{u}}}^{n+1}-{\mathrm{u}}^n}{\Delta t}+[\mathrm{u}\cdot {\nabla }_h\mathrm{u}{]}^{(n+\frac{1}{2})}\right)$$

$${\nabla }_h\cdot {\tilde{\mathrm{u}}}^{n+1}$$

$${\mathrm{F}}^{n+\frac{1}{2}}$$

$$\rho \left(\frac{{\mathrm{u}}^{n+1}-{\tilde{\mathrm{u}}}^{n+1}}{\Delta t}\right)$$

$$[\mathrm{u}\cdot {\nabla }_h\mathrm{u}{]}^{(n+\frac{1}{2})}=\frac{3}{2}{\mathrm{u}}^n\cdot {\nabla }_h{\mathrm{u}}^n-\frac{1}{2}{\mathrm{u}}^{n-1}\cdot {\nabla }_h{\mathrm{u}}^{n-1}.$$

$$-{\nabla }_h{p}^{n+1}+\mu {\nabla }_h^2{\mathrm{u}}^{n+1}+{\mathcal{S}}_h[{\mathrm{X}}^n]{\mathrm{F}}^{n+1}$$

$${\nabla }_h\cdot {\mathrm{u}}^{n+1}$$

$${\mathcal{J}}_h[{\mathrm{X}}^n]{\mathrm{u}}^{n+1}$$

$$\mathcal{F}(t)=-{\oint }_{S_b(t)}\mathrm{n}\cdot \left[-p\mathrm{I}+\mathrm{T}\right]\mathrm{dS},$$

$$\mathcal{F}(t)=-\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{V(t)}\rho \mathrm{u}\mathrm{dV}+{\oint }_{S_{\text{CV}}(t)}\mathrm{n}\cdot \left[-p\mathrm{I}-(\mathrm{u}-{\mathrm{u}}_S)\rho \mathrm{u}+\mathrm{T}\right]\mathrm{dS}-{\oint }_{S_b(t)}\mathrm{n}\cdot (\mathrm{u}-{\mathrm{u}}_S)\rho \mathrm{u}\mathrm{dS},$$

$$\mathcal{F}(t)=-{\int }_{V_{\text{CV}}(t)}\rho \frac{\partial \mathrm{u}}{\partial t}\mathrm{dV}+\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{V_b(t)}\rho \mathrm{u}\mathrm{dV}+{\oint }_{S_{\text{CV}}(t)}\mathrm{n}\cdot \left[-p\mathrm{I}-\mathrm{u}\rho \mathrm{u}+\mathrm{T}\right]\mathrm{dS}.$$

$$\mathcal{M}(t)=-{\oint }_{S_b(t)}\mathrm{r}\wedge \left(\mathrm{n}\cdot \left[-p\mathrm{I}+\mathrm{T}\right]\right)\mathrm{dS},$$

$$\mathcal{M}(t)=-\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{V(t)}(\mathrm{r}\wedge \rho \mathrm{u})\mathrm{dV}+{\oint }_{S_{\text{CV}}(t)}[\mathrm{r}\wedge (-p\mathrm{n}+\mathrm{n}\cdot \mathrm{T})-\mathrm{n}\cdot (\mathrm{u}-{\mathrm{u}}_S)(\mathrm{r}\wedge \rho \mathrm{u})]\mathrm{dS}-{\oint }_{S_b(t)}\mathrm{n}\cdot (\mathrm{u}-{\mathrm{u}}_S)(\mathrm{r}\wedge \rho \mathrm{u})\mathrm{dS}.$$

$$\mathcal{M}(t)=-{\int }_{V_{\text{CV}}(t)}\rho \mathrm{r}\wedge \frac{\partial \mathrm{u}}{\partial t}\mathrm{dV}+\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{V_b(t)}\rho (\mathrm{r}\wedge \mathrm{u})\mathrm{dV}+{\oint }_{S_{\text{CV}}(t)}[\mathrm{r}\wedge (-p\mathrm{n}+\mathrm{n}\cdot \mathrm{T})-(\mathrm{n}\cdot \mathrm{u})\rho (\mathrm{r}\wedge \mathrm{u})]\mathrm{dS}.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{V_b(t)}\rho \mathrm{U}\mathrm{dV}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{V_b(t)}\rho \mathrm{R}\wedge \mathrm{U}\mathrm{dV}$$

$${\int }_{\Omega }\mathrm{f}(\mathrm{x},t)\mathrm{d}\mathrm{x}$$

$$={\int }_{U\subset \Omega }\left[{\int }_{\Omega }\mathrm{F}(\mathrm{s},t)\delta (\mathrm{x}-\mathrm{X}(\mathrm{s},t))\mathrm{d}\mathrm{x}\right]\mathrm{d}\mathrm{s}$$

$$={\int }_{U\subset \Omega }\mathrm{F}(\mathrm{s},t)\left[{\int }_{\Omega }\delta (\mathrm{x}-\mathrm{X}(\mathrm{s},t))\mathrm{d}\mathrm{x}\right]\mathrm{d}\mathrm{s}$$

$$={\int }_U\mathrm{F}(\mathrm{s},t)\mathrm{d}\mathrm{s}.$$

$$\rho \frac{\partial \mathrm{u}}{\partial t}+\rho \nabla \cdot \text{mbox}(\mathrm{u}\mathrm{u})=\nabla \cdot \text{mbox}(-p\mathrm{I}+\mathrm{T})+\mathrm{f},$$

$${\int }_{V_{\text{CV}}(t)}\rho \frac{\partial \mathrm{u}}{\partial t}\mathrm{dV}+{\int }_{V_{\text{CV}}(t)}\rho \nabla \cdot \text{mbox}(\mathrm{u}\mathrm{u})\mathrm{dV}={\int }_{V_{\text{CV}}(t)}\nabla \cdot \text{mbox}\left[-p\mathrm{I}+\mathrm{T}\right]\mathrm{dV}+{\int }_{V_{\text{CV}}(t)}\mathrm{f}\mathrm{dV}.$$

$${\int }_{V_{\text{CV}}(t)}\mathrm{f}\mathrm{dV}={\int }_{V_b(t)}\mathrm{F}\mathrm{dV}={\int }_{V_{\text{CV}}(t)}\rho \frac{\partial \mathrm{u}}{\partial t}\mathrm{dV}-{\int }_{S_{\text{CV}}(t)}\mathrm{n}\cdot \left[-p\mathrm{I}-\mathrm{u}\rho \mathrm{u}+\mathrm{T}\right]\mathrm{dS}.$$

$$y(x,t)=0.125\frac{x+0.03125}{1.03125}\sin \left[2\pi (x-t/T)\right].$$

$$M_{\text{b}}{\mathrm{U}}_{\text{r}}^{n+1}$$

$${\mathrm{I}}_{\text{b}}{\mathrm{W}}_{\text{r}}^{n+1}$$

$$\mathcal{I}(t)=-{\int }_{V_{\text{CV}}(t)}\rho \frac{\partial \mathrm{u}}{\partial t}\mathrm{dV}+{\oint }_{S_{\text{CV}}(t)}\mathrm{n}\cdot \left[-p\mathrm{I}-\mathrm{u}\rho \mathrm{u}+\mathrm{T}\right]\mathrm{dS}=0.$$

$${\mathrm{F}}_{ij}^P=\{\begin{array}{ll}0,& \Vert {\mathrm{x}}_i-{\mathrm{x}}_j\Vert >R_i+R_j+\zeta \text{ or }i=j\\ \frac{c_{ij}}{{ϵ}_P}{\left(\frac{\Vert {\mathrm{x}}_i-{\mathrm{x}}_j\Vert -R_i-R_j-\zeta }{\zeta }\right)}^2\left(\frac{{\mathrm{x}}_i-{\mathrm{x}}_j}{\Vert {\mathrm{x}}_i-{\mathrm{x}}_j\Vert }\right),& \Vert {\mathrm{x}}_i-{\mathrm{x}}_j\Vert \le R_i+R_j+\zeta \end{array}$$

$${\mathcal{F}}_i=\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{V_b^i(t)}\rho {\mathrm{U}}_i\mathrm{dV}-{\int }_{V_b^i(t)}({\mathrm{F}}_i+{\mathrm{F}}_i^g+{\mathrm{F}}_{ij}^P)\mathrm{dV},$$

$$-\nabla p(\mathrm{x})+\mu {\nabla }^2\mathrm{u}(\mathrm{x})+\mathrm{f}(\mathrm{x})=0.$$

$$\mathcal{F}$$

$$\mathcal{M}$$

$$\mathcal{F}$$

Full text

1 Introduction

Fluid-structure interaction (FSI) problems involving moving bodies is a challenging area in the computational fluid dynamics field that has vested the interest of researchers for several decades. FSI modeling has traditionally been carried out in two ways: the body-fitted mesh approach using unstructured grids [1, 2] and the Cartesian grid approach based on the fictitious domain method [3, 4]. Although the body-fitted mesh approach to FSI resolves the fluid-structure interface sharply, it requires complex mesh management infrastructure along with high computational costs for solving linear equations. The fictitious domain method on the other hand extends the fluid equations inside the structure along with some additional non-zero body forcing term. As a result, regular Cartesian grids and fast linear solvers such as fast Fourier transform can be used to solve the common momentum equation of the continua. However, fictitious domain methods tend to smear the fluid-structure interface and hence reduce the solution accuracy near the interface.

One widely used fictitious domain approach to FSI is the immersed boundary (IB) method [5] which was originally proposed by Peskin in the context of cardiac flows [6]. The main advantage to IB methods is that they do not require a body-fitted mesh to model the structure. The immersed body is allowed to freely cut the background Cartesian mesh, making the IB method easy to implement within an existing incompressible flow solver. A Lagrangian force density is computed on structure nodes, which is then transferred to the background grid via regularized delta functions. The use of regularized delta functions diffuses the interface and smears it over a number of grid cells proportional to the width of the delta function. This often leads to discontinuous or noisy derivatives of the velocity and pressure field required to compute surface traction.

The original IB method uses an explicit time stepping scheme and fiber elasticity to compute the additional body forcing inside the region occupied by the structure. The method works fairly well for soft elastic structures, but incurs severe time step restrictions if the stiffness of the material is substantially increased to model rigid bodies. Specialized versions of the original IB method have been developed to model rigid and stiff bodies in an efficient manner like the implicit IB method [7, 8, 9, 10, 11], the direct forcing method [12, 13, 14], and the fully constrained IB method [15, 16, 17]. The implicit IB method requires special solvers like algebraic [18] and geometric multigrid [7, 10] to treat the stiff elastic forces implicitly. The direct forcing and fully constrained IB methods that are primarily used for modeling rigid bodies impose rigidity constraint through Lagrange multipliers. The direct forcing method approximates the Lagrange multiplier with suitable penalty term and solves the fluid-structure equations in a fractional time stepping scheme. This is in contrast to fully constrained IB method where the fluid velocity, pressure, and Lagrange multipliers are solved together. Apart from the obvious advantage of imposing the rigidity constraint exactly rather approximately, the fully constrained IB method can be used for Stokes flow where fractional time stepping schemes do not work. However, special preconditioners are required to solve for the Lagrange multipliers exactly [16, 17], which makes solving the system costly for big three dimensional volumetric bodies.

A common feature of IB methods that are based on Peskin’s IB approach [5] is that they do not require rich geometric information like surface elements and normals to compute the Lagrangian force density. Structure node position (and possibly the node connectivity) information suffices to compute Lagrangian forces. However, there are many other versions of the IB method and Cartesian grid based methods that require additional geometric information to sharply resolve the fluid-structure interface. Examples include the immersed finite element method [19, 20, 21, 22], the immersed interface method [23, 24], the ghost-fluid method [25], and the cut-cell embedded boundary method [26, 27, 28]. For such approaches, the availability of surface elements and normals along with sharp interface resolution makes the surface traction computation easier and smooth (for at least smooth problems).

Often times, the net hydrodynamic forces and torques on immersed bodies are desired rather than point-wise traction values. Force/torque balance equations can be used to compute the hydrodynamic force/torque contribution instead of directly integrating surface traction. In the context of “Peskin-like” IB methods, this implies that one can essentially eliminate surface mesh generation as a post-processing step if only net hydrodynamic forces and torques are desired. This has another added advantage of not using noisy derivatives of velocity and pressure at the interface for evaluating hydrodynamic forces and torques.

In this work we analyze and compare two approaches to computing net hydrodynamic forces and torques on an immersed body. In the first approach we use the Reynolds transport theorem (RTT) to convert the traction integral over an irregular body surface to a traction integral over a regular and simple Cartesian box (that is aligned with grid faces). The RTT is proposed for a moving control volume that translates with an arbitrary velocity to enclose the immersed body at all times. We refer to this approach as the moving CV approach. In the context of locally refined grids, the moving control volume can span a hierarchy of grid levels. In the second approach, hydrodynamic forces and torques are computed using inertia and Lagrange multipliers (approximate or exact) defined in the body region. We refer to this approach as the LM approach and has been used before in [13]. For IB methods, we show that both approaches are equivalent. This is due to a special property of Peskin’s delta functions that makes Lagrangian and Eulerian force density equivalent [5]. We show that both these approaches give smooth forces and suppress spurious oscillations that arise by directly integrating spatial pressure and velocity gradients over the immersed body as reported in the literature [29].

Although application of the RTT on a stationary control volume to evaluate hydrodynamic force is a well known result [30, 31]; its extension to moving control volumes was first proposed by Flavio Noca in 1997 [32]. They were motivated by the task of evaluating net hydrodynamic force on a moving bluff body using DPIV data from experiments. Flavio has also proposed force expressions that eliminate the pressure variable; a quantity not available in DPIV experiments [32, 33]. We do not analyze such expressions in this work, however. In the context of IB method, Bergmann and co-workers [34, 35] have used Flavio’s moving control volume force expressions (involving both velocity and pressure) to compute hydrodynamic forces. They observed spurious force oscillations with a moving control volume approach [34]. We show that by manipulating time derivatives in the original expressions, one can eliminate such spurious oscillations. We also present strategies to mitigate jumps in velocity derivatives in an AMR framework. Such jumps arise when the Cartesian grid hierarchy is regridded and velocity in the new grid hierarchy is reconstructed from the old hierarchy [36]. Lai and Peskin [37] used stationary control volume analysis on a uniform grid to compute steady state hydrodynamic forces on a stationary cylinder for Reynolds numbers between $100-200$ using the IB method. They did not consider time derivative terms in their analysis and temporal hydrodynamic force profiles were reported at steady state. In this work we include time derivative terms for finite Reynolds number flows (but not for steady Stokes flow) in our (moving) control volume analysis.

If point-wise traction values are desired for IB-like methods, there are several recommendations proposed in the literature for smoothing them. Here we list a few of them. Verma et al. [38] have recommended using a “lifted” surface: a surface two grid cell distance away from the actual interface to avoid choppy velocity gradients. This recommendation, based upon their empirical tests using a Brinkman penalization method, can change depending on the smoothness of the problem and the discrete delta function used in IB methods. Goza et al. [39] obtain smooth point-wise force measurements by using a force filtering post-processing step that penalizes inaccurate high frequency stress components. Martins et al. [40] enforce a continuity constraint in the velocity interpolation stencil to reconstruct a second-order velocity field at IB surface points that is discretely divergence-free. They also impose a normal gradient constraint in the pressure interpolation stencil to reconstruct a second-order accurate pressure field at IB points. They have successfully eliminated spurious force oscillations using constrained least-squares stencils. The idea of unconstrained moving least-squares velocity interpolation and force spreading in the context of direct forcing IB method was first proposed by Vanella and Balaras [41]. Lee et al. [29] attribute sources of spurious force oscillations to spatial pressure discontinuities across the fluid-solid interface and temporal velocity discontinuities for moving bodies. They recommend using fine grid resolutions to alleviate spurious oscillations. Their analyses and tests [29] show that grid spacing has a more pronounced effect on spurious force oscillations than computational time step size.

For a range of test cases varying from free-swimming at high Reynolds number to steady Stokes flow, we show that both LM and moving CV methods are in excellent agreement and are able to suppress spurious force oscillations. They do not require any additional treatment such as least-squares stencil (velocity and pressure) interpolation or force filtering beyond simple integration of force balance laws. For moderate to high Reynolds number test cases, we use a direct forcing IB method to estimate Lagrange multipliers, and for Stokes flow we use a fully constrained IB method to compute Lagrange multipliers exactly.

2 Equations of motion

2.1 Immersed boundary method

The immersed boundary (IB) formulation uses an Eulerian description for the momentum equation and divergence-free condition for both the fluid and the structure. A Lagrangian description is employed for the structural position and forces. Let $\bm{\mathrm{x}}=(x_{1},\ldots,x_{d})\in\Omega$ denote fixed Cartesian coordinates, in which $\Omega\subset\mathbb{R}^{d}$ is the fixed domain occupied by the entire fluid-structure system in $d$ spatial dimensions. Let $\bm{\mathrm{s}}=(s_{1},\ldots s_{d})\in U$ denote the fixed material coordinate system attached to the structure, in which $U\subset\mathbb{R}^{d}$ is the Lagrangian curvilinear coordinate domain. The position of the immersed structure occupying a volumetric region $V_{b}(t)\subset\Omega$ at time $t$ is denoted by $\bm{\mathrm{X}}(\bm{\mathrm{s}},t)$. We consider only neutrally buoyant bodies to simplify the implementation; this assumption implies that the fluid and structure share the same uniform mass density $\rho$. The deviatoric stress tensor of fluid, characterized by dynamic viscosity $\mu$, is extended inside the structure to make the momentum equation of both media appear similar. The combined equations of motion for fluid-structure system are [5]

[TABLE]

Eqs. (1) and (2) are the incompressible Navier-Stokes equations written in Eulerian form, in which $\bm{\mathrm{u}}(\bm{\mathrm{x}},t)$ is the velocity, $p(\bm{\mathrm{x}},t)$ is the pressure, and $\bm{\mathrm{f}}(\bm{\mathrm{x}},t)$ is the Eulerian force density, which is non-zero only in the structure region. Interactions between Lagrangian and Eulerian quantities in Eqs. (3) and (4) are mediated by integral equations with Dirac delta function kernels, in which the $d$-dimensional delta function is $\delta(\bm{\mathrm{x}})=\Pi_{i=1}^{d}\delta(x_{i})$. Eq. (3) converts the Lagrangian force density $\bm{\mathrm{F}}(\bm{\mathrm{s}},t)$ into an equivalent Eulerian density $\bm{\mathrm{f}}(\bm{\mathrm{x}},t)$. In the IB literature, the discretized version of this operation is called force spreading. Using short-hand notation, we denote force spreading operation by $\bm{\mathrm{f}}=\bm{\mathrm{\mathcal{S}}}[\bm{\mathrm{X}}]\,\bm{\mathrm{F}}$, in which $\bm{\mathrm{\mathcal{S}}}[\bm{\mathrm{X}}]$ is the force-spreading operator. Eq. (4) determines the physical velocity of each Lagrangian material point from the Eulerian velocity field, so that the immersed structure moves according to the local value of the velocity field $\bm{\mathrm{u}}(\bm{\mathrm{x}},t)$ (Eq. (5)). This velocity interpolation operation is expressed as $\frac{\partial\bm{\mathrm{X}}}{\partial t}=\bm{\mathrm{U}}=\bm{\mathrm{\mathcal{J}}}[\bm{\mathrm{X}}]\,\bm{\mathrm{u}}$, in which $\bm{\mathrm{\mathcal{J}}}[\bm{\mathrm{X}}]$ is the velocity-interpolation operator. It can be shown that if $\bm{\mathrm{\mathcal{S}}}$ and $\bm{\mathrm{\mathcal{J}}}$ are taken to be adjoint operators, i.e. $\bm{\mathrm{\mathcal{S}}}=\bm{\mathrm{\mathcal{J}}}^{*}$, then Lagrangian-Eulerian coupling conserves energy [5].

2.2 Discrete equations of motion

We employ a staggered grid discretization for the momentum and continuity equations (see Fig. 2). More specifically, Eulerian velocity and force variables are defined at face centers while the pressure variable is defined at cell centers. Second-order finite difference stencils are used to spatially discretize the Eulerian equations on locally refined grids [13, 36]. The spatial discretization of various operators are denoted with $h$ subscripts. To discretize equations in time, we take $\Delta t$ as the time step size, and $n$ as the time step number. We use the direct forcing method of Bhalla et al. [13] for moderate to high Reynolds number cases and the fully constrained IB method of Kallemov et al. [16] for Stokes flow cases. The two methods differ in how the Lagrangian force density or Lagrange multipliers are computed. The time integrators are also different for the two methods. Here we briefly describe the discretized equations for both methods. We refer readers to [13, 16, 17] for more details.

For the direct forcing method, the time stepping scheme reads as [13]

[TABLE]

Succinctly, we first solve for a velocity field $\tilde{\bm{\mathrm{u}}}^{n+1}$ and a pressure field $p^{n+\frac{1}{2}}$ as a coupled system by solving Eqs. (6) and (7) simultaneously. The velocity $\tilde{\bm{\mathrm{u}}}^{n+1}$, which is correct in the fluid region but not in the structure region $V_{b}(t)$, is then corrected by estimating the Lagrange multiplier $\bm{\mathrm{F}}^{n+\frac{1}{2}}$ via Eq. (8). Here $\bm{\mathrm{U}}_{b}^{n+1}$ is the desired rigid body velocity of the Lagrangian nodes, and $\bm{\mathrm{X}}^{n+\frac{1}{2}}$ is the midstep estimate of Lagrangian node position. Finally, the Lagrange multiplier $\bm{\mathrm{F}}^{n+\frac{1}{2}}$ is spread on the background grid to correct the momentum in the structure region to $\bm{\mathrm{u}}^{n+1}$. We use Adams-Bashforth to approximate the midstep value of nonlinear convection term $\bm{\mathrm{u}}\cdot\nabla\bm{\mathrm{u}}$ via

[TABLE]

For the fully constrained method, we simultaneously solve for the updated Eulerian velocity $\bm{\mathrm{u}}^{n+1}$ and pressure $p^{n+1}$ at time $t^{n+1}$ along with the Lagrange multiplier $\bm{\mathrm{F}}^{n+1}$. The time stepping scheme reads as [16, 17]

[TABLE]

The coupled system of Eqs. (11)-(13) is solved simultaneously using a preconditioned FGMRES [42] solver. For both methods, we update the Lagrangian node positions $\bm{\mathrm{X}}^{n+1}$ using rigid body translation and rotation. We use Peskin’s 4-point regularized delta functions for the $\bm{\mathrm{\mathcal{S}}}$ and $\bm{\mathrm{\mathcal{J}}}$ operators in all our numerical experiments, unless stated otherwise.

Our finite Reynolds number fluid solver has support for adaptive mesh refinement, and some cases presented in Sec. 5 make use of multiple grid levels (also known as a grid hierarchy). A grid with $\ell$ refinement levels with grid spacing $\Delta x_{0}$, $\Delta y_{0}$, and $\Delta z_{0}$ on the coarsest grid level has minimum grid spacing $\Delta x_{\textrm{min}}=\Delta x_{0}/n_{\text{ref}}^{\ell-1}$, $\Delta y_{\textrm{min}}=\Delta y_{0}/n_{\text{ref}}^{\ell-1}$, and $\Delta z_{\textrm{min}}=\Delta z_{0}/n_{\text{ref}}^{\ell-1}$ on the finest grid level. Here, $n_{\text{ref}}\in\mathbb{N}$ is the refinement ratio. In the present work, the refinement ratio in taken to be the same in each direction, although this is not a limitation of the numerical method. The immersed structure is always placed on the finest grid level. For all of the cases considered in this work a constant time step size $\Delta t=\min(\Delta t^{\ell})$ is chosen, in which the time step size $\Delta t^{\ell}$ on grid level $\ell$ satisfies the convective CFL condition $\Delta t^{\ell}\leq C\min\left(\frac{\Delta x}{\|u_{x}\|_{\infty}},\frac{\Delta y}{\|u_{y}\|_{\infty}},\frac{\Delta z}{\|u_{z}\|_{\infty}}\right)^{\ell}$. In this work, the convective CFL number is set to $C=0.3$ unless otherwise stated.

3 Hydrodynamic force and torque

3.1 Moving control volume method

3.1.1 Hydrodynamic force

Letting $\bm{\mathrm{T}}=\mu\left(\nabla\bm{\mathrm{u}}+\nabla\bm{\mathrm{u}}^{T}\right)$ denote the viscous stress tensor, the net hydrodynamic force is defined to be the force of the fluid on the body:

[TABLE]

in which the integral is taken over the surface of the body $S_{b}(t)=\partial V_{b}(t)$, and $\bm{\mathrm{n}}$ is the unit outward normal to the surface. In practice, evaluating Eq. (14) is inconvenient in numerical experiments because it is often difficult to obtain accurate surface velocity gradients and pressure values. Moreover, evaluating Eq. (14) also requires computational geometry to obtain surface normals and area. Instead we use a control volume approach to compute $\bm{\mathrm{\mathcal{F}}}(t)$ which avoids these requirements.

Let an arbitrary (possibly time dependent) domain $V_{\text{CV}}(t)$ completely surround $V_{b}(t)$, i.e. $V_{b}(t)\subset V_{\text{CV}}(t)$, as shown in Fig. 1. By considering the change in momentum within the control volume $V_{\text{CV}}(t)$ and the net momentum flux at its surface $S_{\text{CV}}(t)=\partial V_{\text{CV}}(t)$, a general expression for the hydrodynamic force on the body can be obtained

[TABLE]

in which $V(t)=V_{\text{CV}}(t)\setminus V_{b}(t)$ is the volume outside the immersed body but inside the CV, and $\bm{\mathrm{u}}_{S}$ is the velocity of the surface over which the (surface) integral is evaluated. In general $\bm{\mathrm{u}}_{S}\neq\bm{\mathrm{u}}$, and $\bm{\mathrm{u}}_{S}$ can be arbitrarily chosen so that the moving CV always encloses the immersed body. The unit normal $\bm{\mathrm{n}}$ in Eq. (15) points outward on $S_{\text{CV}}(t)$, and into $S_{b}(t)$. The above force equation was first derived by Noca [32] and has been used in various experimental [33, 43, 44, 45] and numerical [34, 35, 46, 47] studies. For Cartesian grid based methods, the CV can be chosen as a simple rectangular domain, for which the unit normals on $S_{\text{CV}}(t)$ are aligned with the Cartesian axes. Finally, the integral over $S_{b}(t)$ vanishes for many applications where no-slip ($\bm{\mathrm{u}}_{S}=\bm{\mathrm{u}}$) boundary condition can be chosen for $\bm{\mathrm{u}}_{S}$. Henceforth, we will analyze cases for no-slip boundary conditions.

When Eq. (15) is discretized, the first term requires a discrete approximation of the integral at two separate time instances or at two different locations of the moving control volume. Bergmann et al.  [34] observed spurious force oscillations as a result of this time derivative term. We show that by manipulating the first term using the Reynolds transport theorem, an expression for hydrodynamic force can be obtained that does not require contributions from control volumes at two different spatial locations. The modified equation reads as

[TABLE]

A detailed derivation of Eq. (16) is provided in A. Note that although Eqs. (15) and (16) are equivalent formulas to obtain the hydrodynamic force on an immersed body, their physical interpretations are different. In Eq. (15), the control volume is moving with some prescribed velocity, usually chosen to follow the structure to ensure it is contained within the CV at all times. Hence, one must keep track of the velocity $\bm{\mathrm{u}}_{S}$ of the control surface $S_{\text{CV}}(t)$. Moreover since the time derivative appears outside of the integral over $V(t)$, it requires a discrete approximation of momentum on two time-lagged control volumes. On the other hand in Eq. (16), all integrals over the control volume are evaluated at a single time instance and therefore a discrete evaluation on two separate CVs is never needed. Discretely, the control volume is placed at a new location at time step (${n+1}$) and no information from its previous location at time step $n$ is used. Hence, $\bm{\mathrm{u}}_{S}$ never appears in the calculation. This has the numerical benefit of suppressing force oscillations, which will be shown in Sec. 5.

3.1.2 Hydrodynamic torque

The net hydrodynamic torque on an immersed body is defined to be the net moment of hydrodynamic force exerted by fluid on the body about a given reference point:

[TABLE]

in which $\bm{\mathrm{r}}=\bm{\mathrm{x}}-\bm{\mathrm{x}}_{0}$. The torque is computed with respect to some reference point $\bm{\mathrm{x}}_{0}$, which can be fixed at a location or move with time (e.g. the center of mass of a swimmer). Following Noca’s derivation for the force expression Eq. (15), one can measure the change in angular momentum within a moving control volume to obtain an expression for torque which reads as

[TABLE]

Note that the torque expression in Eq. (18) is slightly different from Eq. 15b given in Bergmann and Iolla [34]. Once again by applying the Reynolds transport theorem to the first term on left-hand side of Eq. (18), we obtain a torque expression involving a control volume contribution at a single spatial location and without any $\bm{\mathrm{u}}_{S}$ terms

[TABLE]

For the derivation of torque expressions (18) and (19), see B.

3.1.3 Numerical integration

We use Riemann summation to evaluate the various integrals in Eqs. (16) and (19) over a moving rectangular control volume. Fig. 2 shows the rectangular control volume marked by its lower and upper coordinates $(x_{L},y_{L})$ and $(x_{U},y_{U})$, respectively.

The arbitrary surface velocity $\bm{\mathrm{u}}_{S}$ is chosen such that the moving CV is forced to align with Cartesian grid faces. This greatly simplifies the evaluation of various terms inside the force and torque integrals. The linear and angular momentum integrals over $V_{b}(t)$ are evaluated in the Lagrangian frame, whereas the rest are computed in the Eulerian frame. The details of these computations are given in C.

3.2 Lagrange multiplier method

The hydrodynamic force and torque on an immersed body can also be computed in an extrinsic manner from force/torque balance laws. Specifically, for an immersed body occupying volume $V_{b}(t)$, the force and torque balance laws read as

[TABLE]

in which $\bm{\mathrm{R}}=\bm{\mathrm{X}}-\bm{\mathrm{X}}_{0}$ is the radius vector from reference point $\bm{\mathrm{X}}_{0}$ to Lagrangian node position $\bm{\mathrm{X}}$, and $\bm{\mathrm{F}}$ is the Lagrange multiplier imposing rigidity constraint as defined in Eq. (3). Note that for the IB method, the net hydrodynamic force $\bm{\mathrm{\mathcal{F}}}$ and torque $\bm{\mathrm{\mathcal{M}}}$ on the body can be readily evaluated in the Lagrangian frame as a part of the solution process without computing any extra terms. We have used this approach in a previous work [13]. We will compare results obtained from both moving CV and LM method in Sec. 5.

3.3 Equivalence of the two methods

Although Eqs. (16) and (20) look different, we now show that the use of (regularized) delta functions in the IB method make them equivalent expressions. To prove this, consider a single body in the domain $\Omega$, occupying a region of space $V_{b}(t)\subset\Omega$. Let $\bm{\mathrm{F}}(\bm{\mathrm{s}},t)$ denote the Lagrange multiplier field defined on Lagrangian nodes, and $\bm{\mathrm{f}}(\bm{\mathrm{x}},t)$ denote its Eulerian counterpart. Here, $\bm{\mathrm{s}}\in U$ and $U\subset\mathbb{R}^{d}$. With continuous and Peskin’s discrete delta functions, the following identity holds

[TABLE]

The above expression is the equivalence of Lagrangian and Eulerian force densities and is a well known result [5, 36]. Starting with the momentum Eq. (1) in Eulerian form

[TABLE]

and integrating it over $V_{\text{CV}}(t)$, we obtain

[TABLE]

Applying the divergence theorem to terms with $\nabla\cdot\mbox{}$, and using Eq. (22) we obtain

[TABLE]

Using the above expression for $\int_{V_{b}(t)}\bm{\mathrm{F}}\,\mathrm{dV}$ with Eq. (20) we get force expression (16). Similarly, one can prove the equivalence of torque expressions (19) and (21) by noting that $\bm{\mathrm{r}}\wedge\nabla\cdot\mbox{}\bm{\mathrm{S}}=\nabla\cdot\mbox{}(\bm{\mathrm{r}}\wedge\bm{\mathrm{S}})$ for any symmetric tensor satisfying $\bm{\mathrm{S}}=\bm{\mathrm{S}}^{T}$.

When multiple bodies exist in the domain, however, there is a subtle difference in the LM and moving CV expressions. The LM expressions in the form of Eqs. (20) and (21) are restricted to individual bodies; therefore in the presence of multiple bodies, $\bm{\mathrm{\mathcal{F}}}$ and $\bm{\mathrm{\mathcal{M}}}$ can be computed separately for each body. For the moving CV method, care must be taken to restrict the control volume to a particular body, i.e. it should not enclose other bodies in its vicinity. Otherwise, $\bm{\mathrm{\mathcal{F}}}$ and $\bm{\mathrm{\mathcal{M}}}$ will contain contribution from multiple bodies. We explore this subtlety with the moving CV method by taking flow past two stationary cylinders in Sec. 5.1.3, and drafting-kissing-tumbling of two sedimenting cylinders in Sec. 5.5. For an AMR framework, the boundary of a control volume can span multiple refinement levels. When velocity is reconstructed from old hierarchy to new, there is generally no guarantee that momentum is conserved because of the inter-level velocity interpolation. This can lead to jumps in $\bm{\mathrm{\mathcal{F}}}$ and $\bm{\mathrm{\mathcal{M}}}$ because of the time derivative terms. Eqs. (20) and (21) are restricted to $V_{b}(t)$, which generally lies on the finest grid level. Therefore, one can expect force and torque calculations to be relatively insensitive to velocity reconstruction operations, which mostly requires intra-level interpolation. We explore this issue with translating plate example in Sec. 5.3.2.

4 Software implementation

We use the IBAMR library [48] to implement the moving control volume method in this work for our numerical tests. IBAMR has built-in support for direct forcing and fully constrained IB methods, among other variants of the IB method. IBAMR relies on SAMRAI [49, 50] for Cartesian grid management and the AMR framework. Solver support in IBAMR is provided by PETSc library [51, 52, 53].

5 Results

5.1 Flow past cylinder

In this section we validate our moving control volume method and Lagrange multiplier method for computing hydrodynamic forces and torques on immersed bodies.

5.1.1 Stationary cylinder

We first consider the flow past a stationary circular cylinder. The cylinder has diameter $D=1$ and is placed in a flow with far-field velocity $\bm{\mathrm{U}}_{\infty}=(U_{\infty},V_{\infty})=(1,0)$. The computational domain is a rectangular channel taken to be of size $18D\times 12D$, with the center of the cylinder placed at $(x,y)=(0,0)$. The domain is discretized by a uniform Cartesian mesh of size $900\times 600$. The aerodynamic drag coefficient $C_{D}=\bm{\mathrm{\mathcal{F}}}\cdot\bm{\mathrm{e}}_{x}/(\rho D\|\bm{\mathrm{U}}_{\infty}\|^{2}/2)$ is calculated numerically in two different ways: via Eq. (16) and via integrating Lagrange multipliers enforcing the rigidity constraint on the cylinder (Eq. (20)). The control volume is taken to be $[-D,1.5D]\times[-D,D]$ and does not move from its initial location. Note that in the case where both the body and the control volume are stationary, Eqs. (15) and (16) equivalent and give the same numerical solution. The density is set to $\rho=1$ and the Reynolds number of the flow is $\text{Re}=\rho U_{\infty}D/\mu=550$. This problem has been studied numerically by Bergmann and Iollo [34] and by Ploumhans and Winckelmans [54]. The temporal behavior of $C_{D}$ matches well with the previous studies [34, 54].

5.1.2 Translating cylinder

Next, we consider the case of a circular cylinder translating with prescribed motion. The parameters used in this case are identical to those of the flow past a stationary cylinder case, except now $\bm{\mathrm{U}}_{\infty}=(U_{\infty},V_{\infty})=(0,0)$. The cylinder is dragged with speed $\bm{\mathrm{U}}_{\text{b}}=(U_{\text{b}},V_{\text{b}})=(-1,0)$. The aerodynamic drag coefficient is computed as $C_{D}=\bm{\mathrm{\mathcal{F}}}\cdot\bm{\mathrm{e}}_{x}/(\rho D\|\bm{\mathrm{U}}_{\text{b}}\|^{2}/2)$. The control volume is initially set to $[-D,1.5D]\times[-D,D]$ and translates to the left every few time steps to ensure that it always contains the cylinder. The density is set to $\rho=1$ and the Reynolds number is $\text{Re}=\rho U_{\text{b}}D/\mu=550$. Periodic boundary conditions are used on all faces of the computational domain. This problem was also studied numerically by Bergmann and Iollo [34].

Fig. 4 shows the time evolution of $C_{D}$ calculated in three different ways, via Eqs. (15) and (16) , and by integrating Lagrange multipliers. The present control volume method matches well with the Lagrange multiplier approach. However, spurious jumps in drag coefficient are seen for the original control volume method outlined by Noca. These spurious oscillations are also present in the computation done by Bergmann and Iollo [34]. We remark that the oscillations seen by Bergmann and Iollo [34] are quantitatively different than the ones presented here in Fig. 4 for comparison, although both use Eq. (15) to compute hydrodynamic force. This can be attributed to differences in the numerical method used to impose constraint: we use a constraint-based immersed boundary method whereas Bergmann and Iollo use a Brinkman penalization method.

5.1.3 Two stationary cylinders

To study the effect of control volume size in the presence of multiple bodies, we consider the case of two stationary circular cylinders each with $D=1$ and placed in a flow with far-field velocity $\bm{\mathrm{U}}_{\infty}=(U_{\infty},V_{\infty})=(1,0)$. The bottom and top cylinders are centered about position $(x,y)=(0,-2D)$ and $(x,y)=(0,2D)$, respectively within a rectangular channel of size $18D\times 12D$. The domain is discretized by a uniform Cartesian mesh of size $900\times 600$ and the Reynolds number is $\text{Re}=\rho U_{\infty}D/\mu=550$ for each cylinder. The density of the fluid is set to $\rho=1$.

Four different (but symmetric) control volume configurations are considered:

Two disjoint CVs located at $[-D,D]\times[-3D,-D]$ and $[-D,D]\times[D,3D]$. 2. 2.

Two CVs located at $[-D,D]\times[-3D,D]$ and $[-D,D]\times[-D,3D]$ that slightly overlap, but do not intersect the other cylinder. 3. 3.

Two CVs located at $[-D,D]\times[-3D,2D]$ and $[-D,D]\times[-2D,3D]$, where each CV holds one full cylinder and half of the second cylinder. 4. 4.

Two CVs located at $[-D,D]\times[-3D,2.7D]$ and $[-D,D]\times[-2.7D,3D]$, where each CV contains both cylinders.

Fig. 5 shows flow visualizations for the four different CV configurations. Fig. 6 shows the drag coefficient over time for each of the four configurations. In the case where the CVs do not overlap (Figs. 5(a) and  6(a)), or when they overlap but do not enclose multiple bodies either partially or fully (Figs. 5(b) and  6(b)) the drag coefficients calculated for the top and bottom cylinders are close to the drag coefficient calculated for a single cylinder considered in Section 5.1.1.

The effect of resizing the control volumes to partially or fully contain other objects is seen in the hydrodynamic drag force measurement. In the case where each CV contains one and a half cylinders (Figs. 5(c) and 6(c)), the drag coefficient deviates significantly from the drag coefficient calculated for a single cylinder. Rather, the computed force is the drag on a combined full and a half cylinder contained within the CV, which is approximately $1.5$ times the drag on a single cylinder. In the case where each CV contains both cylinders (Figs. 5(d) and 6(d)), the measured $C_{D}$ in each control volume is approximately twice the $C_{D}$ measured on a single cylinder. Therefore, in presence of multiple bodies in the domain, care must be taken to restrict the CV to an individual body.

5.2 Oscillating cylinder

In this section we consider cases from previous studies, some of which have reported spurious force oscillations in hydrodynamic drag and lift forces with IB methods. We do not observe such spurious force oscillations using LM and CV methods within an immersed boundary framework.

5.2.1 In-line oscillation

We consider an in-line oscillation of a circular cylinder in a quiescent flow as done in Dütsch et al. [55] and Lee et al. [29]. The cylinder has a diameter $D=1$ and is placed in a domain of size $[-16D,16D]\times[-8D,8D]$ with zero velocity prescribed on all boundaries of the domain. The initial center of mass of the cylinder is placed at $(0,0)$ and its velocity is set to $\bm{\mathrm{U}}_{\text{b}}=\left(-U_{0}\cos(2\pi ft),0\right)$, in which $f$ is the frequency of oscillation. The Reynolds number of the flow is $\text{Re}=\rho U_{0}D/\mu=100$, and the Keulegan-Carpenter number is $\text{KC}=U_{0}/(fD)=5$. The time period of oscillation of the cylinder is given by $T=1/f$. The density of the fluid is set to be $\rho=1$. These parameters are chosen to match those reported in [29, 55].

Three levels of mesh refinement are used, with $n_{\text{ref}}=4$ between each level. The cylinder is embedded in the finest mesh level at all time instances. At the coarsest level, three different mesh sizes are used: $50\times 25$, $100\times 50$, and $200\times 100$, which corresponds to finest grid spacings of $\Delta x_{\textrm{min}}=\Delta y_{\textrm{min}}=0.04D$, $\Delta x_{\textrm{min}}=\Delta y_{\textrm{min}}=0.02D$, and $\Delta x_{\textrm{min}}=\Delta y_{\textrm{min}}=0.01D$, respectively. The computational time step size is chosen to be $\Delta t=0.005D/U_{0}$, matching that of Lee at al.  [29]. This time step size satisfies the convective CFL condition with $C=0.7$, which is found to be stable for all the mesh sizes considered here. A stationary control volume is placed at $[-4D,4D]\times[-2D,2D]$ in order to contain the entire cylinder at all time instances. The drag coefficient is computed as $C_{D}=\bm{\mathrm{\mathcal{F}}}\cdot\bm{\mathrm{e}}_{x}/(0.5\rho U_{0}^{2}D)$.

Fig. 7 shows the time evolution of drag coefficient for the oscillating cylinder. Both CV and LM approaches yield identical results. Even our coarse resolution results are in good agreement with the high resolution results of Lee et al. [29]. Moreover, Lee et al. conducted this case within an immersed boundary framework and observed large spurious force oscillations (see Fig. 14 in [29]) at coarse grid resolutions. They compute drag on the immersed body by evaluating pressure and velocity gradients within the body region $V_{b}(t)$. Only at fine grid resolutions, where the spatial gradients are more accurate, were they able to suppress the spurious force oscillations. Since we do not require spatial gradients of velocity and pressure within the body region $V_{b}(t)$ with our approach, we do not observe such spurious oscillations.

5.2.2 Cross-flow oscillation

Here we consider a cylinder oscillating in the transverse direction with an imposed cross-flow. The cylinder has a diameter $D=1$ and its initial center of mass is placed at $(x,y)=(0,0.2D)$. The cylinder is placed in a domain of size $[-5D,27D]\times[-8D,8D]$ and oscillates in the transverse direction with velocity $\bm{\mathrm{U}}_{\text{b}}=\left(0,-V_{0}\sin(2\pi f_{e}t)\right)$, where $f_{e}$ is the frequency of oscillation. An axial free-stream velocity $U_{\infty}$ is set at the inlet, top, and bottom faces of the computational domain. The transverse traction components are set to zero on the top and bottom boundaries, and the axial and transverse tractions are set to zero at the outflow boundary.

The density of the fluid is set to $\rho=1$ and the Reynolds number based on the free-stream velocity is $\text{Re}=\rho U_{\infty}D/\mu=185$. Letting $f_{0}=0.193U_{\infty}/D$ be the natural shedding frequency for a stationary cylinder, we set $f_{e}/f_{0}=1$. The maximum oscillation velocity of the cylinder is taken to be $0.159V_{0}/(f_{e}D)=0.2$. This case is considered in Lee et al. [29] and Guilmineau and Queutey [56]. A stationary CV is placed at $[-D,D]\times[-2D,2D]$, which contains the cylinder at all time instances (see Fig. 8).

The domain is discretized with three different uniform meshes of sizes $400\times 192$, $800\times 384$, and $1600\times 768$, which corresponds to grid spacings of $\Delta x=\Delta y=0.08D$, $\Delta x=\Delta y=0.04D$, and $\Delta x=\Delta y=0.02D$, respectively. A constant time step size of $\Delta t=0.005D/U_{\infty}$ is used for all the computations, matching that of Lee at al.  [29]. The time step size satisfies the convective CFL condition with $C=0.7$, which is found to be stable for all the mesh sizes considered here. Fig. 9 shows the time evolution of drag coefficient $C_{D}=\bm{\mathrm{\mathcal{F}}}\cdot\bm{\mathrm{e}}_{x}/(0.5\rho U_{\infty}^{2}D)$. Again, both the CV and LM approach produce identical results and do not incur spurious force oscillations even at coarse resolutions as observed by Lee et al. (see Fig. 16 in [29]).

5.2.3 Rotational oscillation

As a last example of this section, we consider the case of a cylinder undergoing a rotational oscillation about its center of mass in a quiescent flow. The diameter of the cylinder is taken to be $D=1$ and is placed in a domain of size $[-20D,20D]\times[-20D,20D]$, with zero velocity prescribed on all computational boundaries. The initial center of mass of the cylinder is placed at $(x,y)=(0,0)$ and it rotates about its center with a velocity $\omega_{b}=A_{m}\sin(2\pi ft)$, in which $f$ is the frequency of oscillation and $T=1/f$ is the time period of oscillation. The density of the fluid is set to be $\rho=1$. The Reynolds number of the flow is $\text{Re}=\rho U_{m}D/\mu=300$, in which $U_{m}=A_{m}D/2$. The cylinder rotates with frequency $f=0.1$ and has maximum angular velocity $A_{m}=10fD$. These parameters are chosen to match Borazjani et al. [57].

The time step size is chosen to be $\Delta t=1\times 10^{-4}T$ and the grid is discretized by a two level mesh, which consists of a coarse mesh of size $512\times 512$, and an embedded fine mesh with refinement ratio $n_{\text{ref}}=4$. The minimum grid spacing at the finest level is $\Delta x_{\textrm{min}}=\Delta y_{\textrm{min}}=0.0195D$. The structure remains on the finest grid level for at all time instances. A stationary control volume is placed at $[-1.01562D,1.01562D]^{2}$ and the torque coefficient is computed as $C_{T}=\bm{\mathrm{\mathcal{M}}}\cdot\bm{\mathrm{e}}_{z}/(0.5\rho U_{m}^{2}D)$. Fig. 10 shows the temporal evolution of $C_{T}$, which is in excellent agreement with sharp-interface CURVIB method results of Borazjani et al. [57].

5.3 Moving plate

In this section we consider the effect of regridding on hydrodynamic force calculations using the control volume approach. In the context of an immersed body AMR framework, regridding occurs when the body has moved some distance, or when new flow features of interest have appeared in the computational domain that require additional mesh refinement to resolve them adequately. We consider a moving plate example to understand the jumps in the drag coefficient due to regridding, and to provide some strategies to mitigate them.

5.3.1 Moving plate on a uniform mesh

First we consider a single level domain case in two spatial dimensions. A finite plate of height $b=1$ is dragged perpendicular to itself with constant velocity $\bm{\mathrm{U}}_{\text{b}}=(U_{\text{b}},V_{\text{b}})=(-1,0)$ in an infinite fluid at rest. The plate is modeled as a thin line of points, separated by grid cell size distance in the transverse direction. The physical domain is a periodic box of dimension $32b\times 22b$ and the domain is discretized by a uniform Cartesian grid of size $1024\times 1024$. The initial location of the control volume is $[-2b,2b]\times[-b,b]$, and it moves with an arbitrary speed to enclose the plate at all time instances. The density of the fluid is set to be $\rho=1$. The Reynolds number of the flow is $\text{Re}=\rho U_{\text{b}}b/\mu=20$. The drag coefficient is calculated as $C_{D}=\bm{\mathrm{\mathcal{F}}}\cdot\bm{\mathrm{e}}_{x}/(\rho U_{\text{b}}^{2}b/2)$. An asymptotic solution $C_{D}\approx 2.09$ was derived by Dennis et al. [58], and this problem was also studied numerically by Bhalla et al. [13].

Fig. 11 shows the time evolution of drag coefficient for the moving plate. We see that the numerical solution obtained by the moving CV computation matches well with the asymptotic value derived in [58] and does not contains any spurious force oscillations or jumps.

5.3.2 Moving plate on an adaptive mesh

Next, we consider the same moving plate example but with locally refined grids. The domain is discretized by a coarse grid of size $256\times 256$. A more refined mesh immediately surrounds the plate with refinement ratio $n_{\text{ref}}=4$, giving the finest level an equivalent grid size of $1024\times 1024$. The mesh is adaptive in the sense that it selectively refines in areas with large velocity gradients and where the immersed body is located in the domain. Apart from the locally refined grids, we use the same parameters of Sec. 5.3.1 for this case.

Three different control volumes are used. First, the CV is set initially to $[-2b,2b]\times[-1.03125b,1.03125b]$ and it translates to the right along with the plate. In this configuration, the CV spans multiple levels of the grid hierarchy. Fig. 12 shows the measured drag coefficient over time. Again, $C_{D}$ evolves towards the asymptotic value derived in [58]. However, there are small jumps in $C_{D}$ over time, which are absent from the uniform mesh case of Sec. 5.3.1. These jumps can be attributed to the mesh hierarchy regridding to follow the moving plate or due to the moving CV itself. To rule out the possibility of jumps due to the motion of control volume, we consider a second CV configuration in which the CV is held stationary all times. The stationary CV again spans multiple grid levels and is big enough to contain the moving plate at all time instances. Fig. 13 shows the CV configuration and measured drag coefficient over time for this CV $=[-2b,-1.03125b]\times[20b,1.03125b]$. The jumps in $C_{D}$ are again observed. In our third configuration, we limit the moving CV on the finest grid level. The CV configuration and $C_{D}$ temporal profile is shown in Fig. 14. The initial location of this control volume is $[-b,0.625b]\times[-1.03125b,1.03125b]$.

From Fig. 14, we note that the jumps due to regridding can be substantially mitigated if the moving CV is restricted to the finest mesh level. The CV translates to the right along with the plate, but remains on the finest mesh level throughout the simulation. There are no longer jumps in the computed drag values for $t\leq 10$.

The jumps in the hydrodynamic forces due to regridding occur due to velocity reconstruction following a regridding operation. A common velocity reconstruction strategy employed in an AMR framework is to first interpolate velocities from the coarser level to the new fine level, and then directly copy velocities from the old fine level in the spatial regions where old and new fine levels intersect. We refer readers to Griffith et al. [36] for details. When we restrict the moving CV to the finest grid level, we reduce errors in momentum change due to velocity reconstruction from coarse to fine level interpolation. The contribution of $\bm{\mathrm{u}}^{n}$ in the time derivative term in Eqs. (16) and (19) can be evaluated either before or after regridding. In our empirical tests, we have observed that evaluating the time derivative term after regridding helps in mitigating the jumps even further (comparison data not shown). In all our results shown above, we evaluate the contribution from $\bm{\mathrm{u}}^{n}$ after regridding, i.e, using the old velocity at the new hierarchy configuration.

5.4 Swimming eel

In this section we demonstrate that the moving control volume approach can be used to determine the hydrodynamic forces and torques on a free-swimming body. We consider a two-dimensional undulating eel geometry, which is adapted from [2, 13]. The eel’s reference frame is aligned with the $x$-axis and in this refrence frame the lateral displacement along $0\leq x\leq L$ over its projected length $L$ is given by

[TABLE]

A backwards-traveling wave of the above form having a time period $T$ causes the eel to self-propel. The swimmer is taken to have a projected length $L=1$, and time period $T=1$. The Reynolds number based on $V_{\textrm{max}}=0.785L/T$, the maximum undulation velocity at the tail tip, is $\text{Re}=\rho V_{\textrm{max}}L/\mu=5609$. The undulations travel in the positive $x$-direction, thereby propelling the eel in the negative $x$-direction. The density of the fluid is taken to be $\rho=1$.

The eel’s total velocity $\bm{\mathrm{U}}_{\text{b}}=\bm{\mathrm{U}}_{\text{r}}+\bm{\mathrm{W}}_{\text{r}}\wedge\bm{\mathrm{R}}+\bm{\mathrm{U}}_{\text{k}}$ in the Lagrangian frame is given by the following components: its rigid linear $\bm{\mathrm{U}}_{\text{r}}$ and angular $\bm{\mathrm{W}}_{\text{r}}$ center of mass velocities and its deformational velocity $\bm{\mathrm{U}}_{\text{k}}$, which we assume to have zero net linear and angular momentum. The self-propulsion velocities are obtained using conservation of linear and angular momentum in the body domain $V_{b}(t)$

[TABLE]

in which $M_{\text{b}}$ and $\bm{\mathrm{I}}_{\text{b}}$ are the mass and moment of inertia tensor of the body, respectively. Having obtained these rigid body velocities, the body velocity at time $t^{n+1}$ required in Eq. (8) is obtained as $\bm{\mathrm{U}}_{\text{b}}^{n+1}=\bm{\mathrm{U}}_{\text{r}}^{n+1}+\bm{\mathrm{W}}_{\text{r}}^{n+1}\wedge\bm{\mathrm{R}}^{n+\frac{1}{2}}+\bm{\mathrm{U}}_{\text{k}}^{n+1}$, in which $\bm{\mathrm{R}}^{n+\frac{1}{2}}=\bm{\mathrm{X}}^{n+\frac{1}{2}}-\bm{\mathrm{X}}^{n+\frac{1}{2}}_{0}$ is the radius vector from (an estimated) midstep center of mass $\bm{\mathrm{X}}_{0}^{n+\frac{1}{2}}$ to (an estimated) midstep Lagrangian node position $\bm{\mathrm{X}}^{n+\frac{1}{2}}$. We refer readers to [13] for more details.

The fully periodic domain is taken to be of size $8L\times 4L$ and is discretized with a three-level hierarchy of Cartesian grids. The size of the coarsest grid is $128\times 64$ grid cells and $n_{\text{ref}}=4$ is taken for subsequent finer grids. Hence, the finest grid, with spacing equivalent to that of a uniform mesh of size $2048\times 1024$, embeds the undulatory swimmer at all times. A time step size of $\Delta t=1\times 10^{-4}\;T$ is employed. The head of the swimmer is initially centered at $(x,y)=(0,0)$ and its body extends in the positive $x$-direction. The CV is initially located at $[-1.02L,1.0425L]\times[-0.7075L,0.73L]$, which encompasses the entire swimmer. The CV is allowed to span multiple grid levels. Whenever the eel’s center of mass translates a distance $\Delta x$, the CV moves with velocity $\bm{\mathrm{u}}_{S}=(-\Delta x/\Delta t,0)$ in order to remain aligned with the grid lines. Fig. 15 shows the vortical structures generated by the eel at four separate time instances, along with the locations of the moving control volume.

Fig. 16 shows the time evolution of axial ($U_{r}=\bm{\mathrm{U}}_{\text{r}}\cdot\bm{\mathrm{e}}_{x}$) and lateral ($V_{r}=\bm{\mathrm{U}}_{\text{r}}\cdot\bm{\mathrm{e}}_{y}$) swimming velocities, along with the rotational velocity $W^{z}_{r}=\bm{\mathrm{W}}_{\text{r}}\cdot\bm{\mathrm{e}}_{z}$. The eel is shown to travel in the $-\bm{\mathrm{e}}_{x}$ direction, eventually reaching a steady state speed. The angular velocity oscillates about a zero mean value, while the lateral velocity has small non-zero mean due to initial transients. Fig. 16 shows the time evolution of net axial ($\mathcal{F}_{x}=\bm{\mathrm{\mathcal{F}}}\cdot\bm{\mathrm{e}}_{x}$) and lateral ($\mathcal{F}_{y}=\bm{\mathrm{\mathcal{F}}}\cdot\bm{\mathrm{e}}_{y}$) forces acting on the eel’s body, along with the net torque $\mathcal{M}_{z}=\bm{\mathrm{\mathcal{M}}}\cdot\bm{\mathrm{e}}_{z}$, which is measured from the eel’s center of mass. The forces and torque are computed using the moving CV approach, although identical estimates are obtained from the LM approach (data not shown). Both $\bm{\mathrm{\mathcal{F}}}$ and $\bm{\mathrm{\mathcal{M}}}$ oscillate about a mean value of zero, which is expected during free-swimming as there are no external forces and torques applied on the swimmer.

For an object initially at rest in a periodic and quiescent fluid, the net linear momentum over the entire computational domain should remain zero [13], i.e., $\bm{\mathrm{\mathcal{P}}}(t)=\int_{\Omega}\rho\bm{\mathrm{u}}\,\mathrm{dV}=\bm{\mathrm{0}}$. Similarly, the net angular momentum of the system should also remain zero at all times, i.e., $\bm{\mathrm{\mathcal{L}}}(t)=\int_{\Omega}\rho\bm{\mathrm{r}}\wedge\bm{\mathrm{u}}\,\mathrm{dV}=\bm{\mathrm{0}}$. In other words, all of the momentum generated due to the eel’s vortex shedding should be redistributed to the eel’s translational and rotational motion. Moreover, the change in linear momentum of the body should be equal to the net force on the body during free-swimming. Hence, the net force on the body should be given by $\bm{\mathrm{\mathcal{F}}}(t)=\frac{{\rm d}}{{\rm d}t}\int_{V_{b}}\rho\bm{\mathrm{u}}\,\mathrm{dV}$, implying that

[TABLE]

The above statement also implies the conservation of linear momentum for the control volume. Eq. (28) also implies that the sum of the Langrange multipliers enforcing the rigidity constraint within $V_{b}(t)$ is zero during free-swimming (see Eq. (24)). Fig. 16 shows the temporal evolution of ${\mathcal{P}}_{x}=\bm{\mathrm{\mathcal{P}}}\cdot\bm{\mathrm{e}}_{x}$, ${\mathcal{P}}_{y}=\bm{\mathrm{\mathcal{P}}}\cdot\bm{\mathrm{e}}_{y}$, $\mathcal{I}_{x}=\bm{\mathrm{\mathcal{I}}}\cdot\bm{\mathrm{e}}_{x}$, and $\mathcal{I}_{y}=\bm{\mathrm{\mathcal{I}}}\cdot\bm{\mathrm{e}}_{y}$; it is indeed seen that these quantities are nearly zero for all time instances. The slight increase in ${\mathcal{P}}_{x}$ is attributed to spatial and temporal discretization errors, whereas the jumps in $\mathcal{I}_{x}$ correspond to time steps at which regridding occurs. Similar observations are made for angular momentum conservation for the entire system (data not shown here).

5.5 Drafting, kissing, and tumbling

In this section we simulate the dynamic interactions between two sedimenting cylindrical particles and use the moving control volume and Lagrange multiplier approaches to determine the hydrodynamic forces. The cylinders are identically shaped with diameter $D=0.2$ cm and are placed in a domain of size $[-5D,5D]\times$ [0, 40D], with zero velocity prescribed on the left and right boundaries, and with axial and transverse tractions set to zero at the top and bottom boundaries. The density and viscosity of the fluid are set to $\rho=1.0$ g/cm3 and $\mu=0.01$ g/(cm $\cdot$ s), respectively. Each particle is subject to a gravitational body force $\bm{\mathrm{F}}^{g}=-(\rho_{s}-\rho)gV_{p}\bm{\mathrm{e}}_{y}$, where $g=980$ cm/s2 is the gravitational constant, $\rho_{s}=1.01\rho$ is the density of the solid, and $V_{p}=\pi(D/2)^{2}$ is the volume of each particle. This is realized through an Eulerian body force $\bm{\mathrm{f}}^{g}$ added to the right-hand side the momentum Eq. (1), which is nonzero only in the particle domains. Similar to the free-swimming eel case, the particles’ translational and rotational velocities are obtained via Eqs. (26) and (27).

The domain is discretized with a two-level hierarchy of Cartesian grids. The size of the coarsest grid is $64\times 256$ grid cells and $n_{\text{ref}}=4$ is taken for subsequent finer grids. Hence, the finest grid, with spacing equivalent to that of a uniform mesh of size $256\times 1024$, embeds each particle at all times. The minimum grid spacing on the finest level is $\Delta x_{\textrm{min}}=\Delta y_{\textrm{min}}=0.0390625D$. A time step size of $\Delta t=5\times 10^{-4}$ s is used. Particle $1$ is placed with initial center of mass $\bm{\mathrm{x}}_{1}=(X_{p},Y_{p})=(-0.005D,36D)$, while particle $2$ is placed below particle $1$ with initial center of mass location at $\bm{\mathrm{x}}_{2}=(X_{p},Y_{p})=(0,34D)$. Under these conditions, the two particles start to accelerate downwards due to gravity. Particle $1$ travels through a low pressure wake created by the leading particle $2$, which causes particle $1$ to fall faster; this stage is called drafting. Eventually, particle $1$ catches up to and nearly contacts particle $2$, a process termed as kissing in literature. This kissing stage is unstable and eventually the particles are left to tumble separately. The parameters here are chosen to match with previous numerical studies on drafting, kissing, and tumbling done by Feng et al. [59], Jafari et al. [60], and Wang et al. [61].

Artificial repulsive forces are added to avoid numerical issues due to overlapping particles. The functional form of this force on particle $i$ due to particle $j$ is given by

[TABLE]

in which $R_{i}=R_{j}=R$ is the radius of both particles, $c_{ij}=\rho\pi R^{2}g$ is a force scale parameter, $\epsilon_{P}=2.0$ g cm/s2 is a stiffness parameters for collisions, and $\zeta=\Delta y_{\textrm{min}}$ is a mesh threshold parameter indicating how far away the two particles need to be in order to feel particle-particle interaction force. This particular repulsive force was used by Feng et al. [59]. No repulsive force between the particle and the wall are used for the case considered here. Similar to the gravitational force, the particle interaction force is realized through an Eulerian body force $\bm{\mathrm{f}}_{P}$ added to the right-hand side of the momentum Eq. (1), which is nonzero only in the particle domains.

The strategy for placing the control volume for this example is different than the previous examples. Rather than setting the CV in motion with some prescribed velocity to enclose the body, the CVs are chosen to surround each cylinder based on its center of mass location: $[X_{p}-5\Delta x_{\textrm{min}},X_{p}+5\Delta x_{\textrm{min}}]\times[Y_{p}-5\Delta y_{\textrm{min}},Y_{p}+5\Delta y_{\textrm{min}}]$. This means that when the particles are close to each other during the kissing stage, each CV contains the second body partially, leading to inaccurate force measurements. This is a limitation of using simple rectangular control volumes. Fig. 17 shows the vortical structures generated by the two particles at four separate time instances, along with the locations of the moving control volumes.

The expression for hydrodynamic force based on the Lagrange multiplier method needs to be modified to account for the presence of additional body forces in the solid region (additional to the Lagrange multiplier constraint forces). These forces need to be added to the right-hand side of Eq. (20), yielding

[TABLE]

in which $\bm{\mathrm{\mathcal{F}}}_{i}$, $\bm{\mathrm{U}}_{i}$, $\bm{\mathrm{F}}_{i}$, $\bm{\mathrm{F}}_{i}^{g}$and $V_{b}^{i}(t)$ are the net hydrodynamic force, center of mass velocity, Lagrange multiplier force, gravitational force, and domain of particle $i$, respectively. The expression for hydrodynamic force based on control volume analysis remains unchanged in the presence of the additional body forces in the solid region, and Eq. (16) remains valid in this scenario111This is because $\int_{V_{\text{CV}}(t)}\bm{\mathrm{f}}\,\mathrm{dV}$ that contains contribution of all body forces in Lagrangian domain is evaluated via the right-hand side of Eq. (24).

Figs. 18 and 18 show the time evolution of the centers of mass $(X_{p},Y_{p})$ of both particles. Particle $1$ gradually approaches particle $2$ up until near $t=2.5$ s when the particles kiss and eventually separate. Upon separation, both particles over time reach a terminal velocity. This temporal behavior matches well with the results of Jafari et al. [60]. Figs. 18 and 18 show the time evolution of net axial ($\mathcal{F}_{x}=\bm{\mathrm{\mathcal{F}}}\cdot\bm{\mathrm{e}}_{x}$) and lateral ($\mathcal{F}_{y}=\bm{\mathrm{\mathcal{F}}}\cdot\bm{\mathrm{e}}_{y}$) forces acting on each particle, calculated using both the control volume and Lagrange multiplier approaches. Between $t=1$ s and $t=2.5$ s, the particles are close to each other and it is not possible to create rectangular CVs that contain only a single particle (see insets in Figs. 17(b) and  17(c)). Hence, the CV force calculations during this time period are inaccurate. However, the forces calculated by the LM method remain accurate at all times. Outside of this time period, both the CV and LM approaches are in excellent agreement. Eventually, the net hydrodynamic force on each particle is $-\bm{\mathrm{F}}^{g}$, indicating a terminal velocity has been achieved.

5.6 Stokes flow

In the previous sections, we considered finite Reynolds number cases simulated using a direct forcing IB method in which Lagrange multipliers were approximated in the body domain. Here we consider a fully constrained IB method in which we compute Lagrange multipliers exactly. We consider Stokes flow examples here, although the fully constrained method also work equally well at finite Reynolds numbers as shown in [16].

For steady Stokes flow in the absence of inertia ($\rho=0$), the momentum equation reads as

[TABLE]

Since Eq. (31) is a steady state problem, it cannot be solved numerically with the split fluid-structure solver described in [13]. Rather, the discretized system of Eqs. (11)-(13) are solved by the monolithic fluid-structure solver described in [16] to obtain a numerical solution to the constrained Stokes system. The rigidity constraint is enforced on the surface of the body, and the body is discretized only by surface nodes and not by a volumetric mesh. This is because enforcing the rigidity constraint on the surface also imposes rigid body motion inside the body for Stokes flow. Setting the inertial terms to zero in Eqs. $\eqref{eq_OurForce}$ and $\eqref{eq_OurTorque}$ yield

[TABLE]

Hence, the net hydrodynamic force and torque on an object in Stokes flow is simply the pressure and viscous fluxes through the control surface. For the LM method, the Lagrange multiplier force density $\bm{\mathrm{F}}(\bm{\mathrm{X}})$ is computed (exactly) on body surface $\bm{\mathrm{X}}\in S_{b}$, and the net force and torque is given by

[TABLE]

5.6.1 Flow between two concentric shells

We first consider the case of two concentric shells, which was studied numerically by Kallemov et al. [16]. The inner and outer shell have geometric radii $R_{1}^{g}=1.807885$ and $R_{2}^{g}=4R_{1}^{g}$, respectively. The computational domain is a cube of size $[0,L]^{3}=[0,4.15R_{2}^{g}]^{3}$, which is discretized by a uniform grid of size $60^{3}$. The center of both shells is placed at $(x,y,z)=(L/2,L/2,L/2)$. Uniform velocity $\bm{\mathrm{u}}=(1,0,0)$ is prescribed on each wall of the computational domain, and the inner and outer shell are set to have rigid–body velocity $(U_{1},V_{1},W_{1})=(0,0,0)$ and $(U_{2},V_{2},W_{2})=(1,0,0)$, respectively. The inner shell is discretized with $42$ surface markers while the outer shell is discretized with $642$ surface markers to ensure that the markers are about $2$ grid cells apart. The viscosity is set to $\mu=1$.

Each spherical shell has an effective hydrodynamic radius $R^{h}$ due to the immersed boundary kernel used to discretize the delta-function [62]. For the 6-point kernel considered by Kallemov et al., it was found that $R_{1}^{h}=1.22R_{1}^{g}$ and $R_{2}^{h}=0.96R_{2}^{g}$ for the numerical parameters chosen here [16]. It was also found that as both the Eulerian and Lagrangian meshes are refined, $R^{h}/R^{g}\rightarrow 1$. These hydrodynamic radii can be used in the analytical expression for the drag on the inner sphere [63], given by

[TABLE]

in which $K=(1-\lambda^{5})/\alpha$, $\alpha=1-9\lambda/4+5\lambda^{3}/2-9\lambda^{5}/4+\lambda^{6}$ and $\lambda=R_{1}^{h}/R_{2}^{h}$.

Depending on the control volume size, the drag on either the inner shell, or on both shells can be obtained. First, a CV of dimension $[L/2-1.659R_{1}^{g},L/2+1.659R_{1}^{g}]^{3}$ is chosen to surround the inner shell, but to exclude the outer shell. Next, a CV of dimension $[L/2-1.383R_{2}^{g},L/2+1.383R_{2}^{g}]$ is chosen to include both inner and outer shells. Fig. 19 shows the configuration of the concentric shells/spheres and the two control volumes.

Table 1 shows the drag measurements from the analytical expression, from integrating surface Lagrange multipliers, and from the control volume analysis. The middle column shows that all three methods are in agreement for the drag on the inner shell. Moreover, the last column shows that the combined drag on both inner and outer shells is the same when computed from Lagrange multipliers and control volume analysis.

5.6.2 Single rotating shell

Next, we consider a single shell rotating in a bounded domain with no exterior flow. The shell is taken to be the same as the inner sphere in the previous example, with geometric radius $R_{1}^{g}=1.807885$ in a computational domain of size $[0,L]^{3}=[0,16.6R_{1}^{g}]^{3}$. There is no outer shell in this example. The center of the shell is placed at the centroid of the cube $(x,y,z)=(L/2,L/2,L/2)$ and $\bm{\mathrm{u}}=(0,0,0)$ is set at all computational boundaries. Periodic boundary conditions were also used and yielded nearly identical torque measurements (data not shown). The viscosity is set to $\mu=1$. The shell rotates about a diameter with angular velocity $\bm{\mathrm{\omega}}=(0,1,0)$. In an unbounded flow at rest, Faxén’s law states that the torque on the sphere by the fluid is

[TABLE]

in which the hydrodynamic radius of the sphere $R_{1}^{h}$ is used [63]. Although the domain in the numerical method is bounded, we still get decent agreement between our numerical results and Eq. (37).

A CV of dimension $[L/2-1.659R_{1}^{g},L/2+1.659R_{1}^{g}]^{3}$ is chosen to surround the shell. Three different grid sizes are used to discretize the domain: $60^{3}$, $120^{3}$, and $240^{3}$, which corresponds to $42$, $162$, and $642$ surface markers on the shell respectively (to ensure that the markers are approximately 2 grid cells apart).

Table 2 shows the torque measurements $\bm{\mathrm{\mathcal{M}}}\cdot\bm{\mathrm{e}}_{y}$ from the analytical expression, from integrating the moments of surface Lagrange multipliers, and from the control volume analysis, for the three different grid resolutions. As expected, the LM and CV measured torque do not match exactly with the analytical expression (presumably because of finite domain effects). However, the LM and CV torque values are in excellent agreement with each other.

Conclusions

In the present study we presented a moving control volume (CV) approach to compute the net hydrodynamic forces and torques on a moving body immersed in a fluid. This approach does not require evaluation of (possibly) discontinuous spatial velocity or pressure gradients within or on the surface of the immersed body. The analytical expressions for forces and torques were modified from those initially presented in [32], and this modification has been shown to eliminate spurious jumps in drag [34]. Our implementation treats the control volume as a rectangular box whose boundary is forced to remain on grid lines, which greatly simplifies the evaluation of surface integrals.

The approach is shown to accurately compute the forces and torques on a wide array of fluid-structure interaction problems, including flow past stationary and moving objects, Stokes flow, and high Reynolds number free-swimming. Spurious momentum gain or loss due to adaptive mesh refinement can produce jumps in the computed forces in the CV approach, although forcing the CV to remain on the finest grid level can ameliorate this issue.

We also show the equivalence between the Lagrange multiplier (LM) approach and the CV approach. The main advantage of the CV approach over the LM approach is that it is applicable to situations where explicit Lagrange multipliers are not available, for example, in the embedded boundary/cut-cell approach to FSI.

The control volume approach implemented here assumes a no-slip boundary condition on the fluid-structure interface. However, a generalization for transpiration boundary conditions can be derived as well (see Eq. (44)). Use of such a boundary condition would required richer geometric information and data structures to evaluate surface quantities of the immersed body. Finally, our approach can be easily extended to cases where additional body forces are present in the momentum equation.

Acknowledgements

A.P.S.B and N.N acknowledge helpful discussions related to software design in IBAMR with Boyce E. Griffith (UNC-Chapel Hill) over the course of this work. N.N, N.A.P, and A.P.S.B acknowledge computational resources provided by Northwestern University’s Quest high performance computing service. N.N acknowledges research support from the National Science Foundation Graduate Research Fellowship Program (NSF award DGE-1324585). N.A.P acknowledges support from the National Science Foundation (NSF award SI2-SSI-1450374). A.P.S.B and H.J acknowledge support from the U.S. Department of Energy, Office of Science, ASCR (award number DE-AC02-05CH11231). A part of this work was carried at UNC-Chapel Hill for which A.P.S.B gratefully acknowledges support from awards NIH HL117163 and NSF ACI 1450327 (awarded to Boyce E. Griffith).

Appendix A Derivation of the new hydrodynamic force expression

Here we present a detailed derivation of Eq. (16) from Eq. (15). The ultimate goal is to obtain an expression for hydrodynamic force which involves integral contributions from a single CV rather than two time-lagged CVs. The Reynolds transport theorem (RTT) [30, 31] gives an expression for the time derivative of an arbitrary quantity $\bm{\mathrm{q}}(\bm{\mathrm{x}},t)$ on a time dependent region $\Omega(t)$

[TABLE]

in which $\bm{\mathrm{u}}_{S}$ is the velocity and $\bm{\mathrm{n}}$ is the outward pointing unit normal vector of the boundary $\partial\Omega(t)$. Applying the RTT to Eq. (15) yields the expression

[TABLE]

Recall that $V(t)=V_{\text{CV}}(t)\setminus V_{b}(t)$, where $V_{\text{CV}}(t)$ is the entire control volume and contains the body domain $V_{b}(t)$; $S_{\text{CV}}(t)=\partial V_{\text{CV}}(t)$ is the boundary of the CV, and $S_{b}(t)=\partial V_{b}(t)$. Hence, $\partial V(t)=S_{\text{CV}}(t)\cup S_{b}(t)$. and the second integral can be split into two boundary integrals

[TABLE]

which can be simplified to obtain

[TABLE]

Recall that in Eq. (41), the unit normal vector points outward on $S_{\text{CV}}(t)$ and inward on $S_{b}(t)$ since the integral is considered with respect to the boundary of $V(t)$. Next, notice that by definition $V_{\text{CV}}(t)=V(t)\cup V_{b}(t)$ is the union of disjoint regions. Hence, the integral over $V(t)$ can be split into $\int_{V(t)}\bm{\mathrm{q}}\,\mathrm{dV}=\left(\int_{V_{\text{CV}}(t)}-\int_{V_{b}(t)}\right)\bm{\mathrm{q}}\,\mathrm{dV}$. Applying the split to Eq. (41) yields

[TABLE]

Finally, we can apply the Reynolds transport theorem (Eq.  (38)) to the integral over $V_{b}(t)$ above. Letting $\bm{\mathrm{N}}$ be the outward pointing unit normal vector to $S_{b}(t)$, we obtain

[TABLE]

Substituting the fact that $\bm{\mathrm{N}}=-\bm{\mathrm{n}}$ yields a general expression for the hydrodynamic force on an immersed body

[TABLE]

Eq. (44) can be further simplified if we assume no-slip boundary conditions at the fluid-structure interface by setting $\bm{\mathrm{u}}_{S}=\bm{\mathrm{u}}$, which gives

[TABLE]

This is the net hydrodynamic force expression as written in Eq. (16).

Appendix B Derivation of the new hydrodynamic torque expression

Conservation of angular momentum for a material volume $V_{m}(t)$ (a volume that moves with the local fluid velocity) that shares its boundary with an arbitrary moving control volume $V(t)$ at time $t$ can be written as [30, 31]

[TABLE]

in which $\bm{\mathrm{r}}=\bm{\mathrm{x}}-\bm{\mathrm{x}}_{0}$, with $\bm{\mathrm{x}}_{0}$ as a reference point for computing torques, and $\Sigma=-p\bm{\mathrm{I}}+\bm{\mathrm{T}}$. Letting $\bm{\mathrm{\mathcal{M}}}(t)$ be the torque exerted by the fluid on the body and noticing that $\oint_{S_{b}(t)}\bm{\mathrm{r}}\wedge[\bm{\mathrm{n}}\cdot\Sigma]\,\mathrm{dS}=-\bm{\mathrm{\mathcal{M}}}(t)$, we have

[TABLE]

Using the RTT, the integral of an arbitrary quantity $\Phi$ over the material volume $V_{m}(t)$ can be related to integral over arbitrary volume $V(t)$ with surface velocity moving with $\bm{\mathrm{u}}_{S}$ as

[TABLE]

Using Eq. (48) with $\Phi=\bm{\mathrm{r}}\wedge\rho\bm{\mathrm{u}}$, the expression for torque becomes

[TABLE]

Finally, by manipulating the term derivative term in Eq. (49) using the RTT we get an expression for torque on an immersed body as

[TABLE]

Appendix C Numerical discretization

Here we describe the discrete evaluation of Eqs. (16) and (19) to obtain the net hydrodynamic force and torque on an immersed body. For notational simplicity, we present the discretized equations in two spatial dimensions. An extension to three spatial dimensions is straightforward. A discrete grid covers the physical domain $\Omega$ with mesh spacing $\Delta x$ and $\Delta y$ in each direction. The position of each grid cell center is given by $\bm{\mathrm{x}}_{i,j}=\left(x_{i,j},y_{i,j}\right)$. For a given cell center, $\bm{\mathrm{x}}_{i-\frac{1}{2},j}$ denotes the physical location of the cell face that is half a grid space away from $\bm{\mathrm{x}}_{i,j}$ in the negative $x$-direction, i.e. $\bm{\mathrm{x}}_{i-\frac{1}{2},j}=\left(x_{i,j}-\frac{\Delta x}{2},y_{i,j}\right)$. Similarly $\bm{\mathrm{x}}_{i,j-\frac{1}{2}}$ denotes the physical location of the cell face that is half a grid cell away from $\bm{\mathrm{x}}_{i,j}$ in the negative $y$-direction, i.e. $\bm{\mathrm{x}}_{i,j-\frac{1}{2}}=\left(x_{i,j},y_{i,j}-\frac{\Delta y}{2}\right)$. The discrete approximations described here are also valid when adaptive mesh refinement is used, although volume weights need to be appropriately modified for different velocity components.

Let $t^{n}$ be the time at time step $n$. After stepping forward from time $t^{n}$ to $t^{n+1}=t^{n}+\Delta t$, a pressure solution is obtained at cell centers $p_{i,j}^{n+1}=p\left(\bm{\mathrm{x}}_{i,j},t^{n+1}\right)$, while velocity components are obtained at cell faces: $u_{i-\frac{1}{2},j}^{n+1}=u\left(\bm{\mathrm{x}}_{i-\frac{1}{2},j},t^{n+1}\right)$ and $v_{i,j-\frac{1}{2}}^{n+1}=v\left(\bm{\mathrm{x}}_{i,j-\frac{1}{2}},t^{n+1}\right)$. These are the only Eulerian quantities needed to evaluate the discrete approximations to Eqs. (16) and (19). Only rectangular control volumes are considered in the present work. A CV is described by its lower left and upper right corners: letting $(x_{L},y_{L})$ and $(x_{U},y_{U})$ denote the lower and upper corners respectively, the control volume is defined to be the Cartesian product of intervals $V_{\text{CV}}(t)=\left\{\bm{\mathrm{x}}\in\Omega\mid\bm{\mathrm{x}}\in[x_{L},x_{U}]\times[y_{L},y_{U}]\right\}$. Moreover, $S_{\text{CV}}(t)$ is forced to remain on grid lines and it not allowed to cross into the interior of grid cells. This greatly simplifies the required numerical approximations. Refer to Fig. 2 for a sketch of the control volume configuration over a staggered mesh discretization. Let the control volume and surface at a time instance $t^{n+1}$ be denoted by $V_{\text{CV}}^{n+1}=V_{\text{CV}}(t^{n+1})$ and $S_{\text{CV}}^{n+1}=S_{\text{CV}}(t^{n+1})$, respectively.

C.1 Discrete approximation to surface integrals

The control surface $S_{\text{CV}}^{n+1}$ is composed of four segments (eight faces in $3$D) denoted by $\mathcal{B}$, ${\mathcal{L}}$, $\mathcal{T}$, and ${\mathcal{R}}$ in Fig. 2. Consequently, computing surface normals on each of these segments is simple, e.g for the bottom segment $\mathcal{B}^{n+1}$, $\bm{\mathrm{n}}=-\bm{\mathrm{e}}_{y}$. The discretized surface integral of a quantity $\bm{\mathrm{\Phi}}$ over $S_{\text{CV}}^{n+1}$ is simply the sum over these four segments

[TABLE]

Moreover, it is sufficient to show the discrete approximation to the surface integral over a single segment since the contribution from the other three segments are computed analogously. Over the bottom surface $\mathcal{B}$, the discretization of each term is given by

[TABLE]

Fig. 20 shows a schematic of the pressure and velocity values required to evaluate Eqs. (52), (53), and (54). Evaluating the surface integral on $\mathcal{B}$ in the torque calculation about a point $\bm{\mathrm{x}}_{0}$ is done by computing $\bm{\mathrm{r}}_{i,j-\frac{1}{2}}=\bm{\mathrm{x}}_{i,j-\frac{1}{2}}-\bm{\mathrm{x}}_{0}$ and evaluating the cross product between $\bm{\mathrm{r}}_{i,j-\frac{1}{2}}$ and the integrand.

C.2 Change in control volume momentum

A discretization of the the time derivative term in Eq. (16) for the integral over $V_{\text{CV}}(t)$ at time step $n+1$ is given by

[TABLE]

where a discrete approximation to the total linear momentum within $V_{\text{CV}}^{n+1}$ can be written as

[TABLE]

Here, $\Delta V=\frac{\Delta x\Delta y}{2}$ when either $(i-\frac{1}{2},j)\in S_{\text{CV}}^{n+1}$ or $(i,j-\frac{1}{2})\in S_{\text{CV}}^{n+1}$, and $\Delta V=\Delta x\Delta y$ otherwise, to ensure that $\sum_{(i-\frac{1}{2},j)}\Delta V_{i-\frac{1}{2},j}=\sum_{(i,j-\frac{1}{2})}\Delta V_{i,j-\frac{1}{2}}=\left|V_{\text{CV}}^{n+1}\right|$, the volume of the CV.

In the original hydrodynamic force formula Eq. (15) introduced by Noca [32], a discretization of the momentum term is given by

[TABLE]

Notice that Eqs. (55) and (57) are nearly identical, although the former only requires an evaluation over a single CV, while the latter requires an evaluation over two time-lagged CVs.

The analogous term in the torque calculation Eq. (19) is discretized differently. Each $u$ velocity location $\bm{\mathrm{x}}_{i-\frac{1}{2},j}$ is looped over and an approximation to $v$ is computed. Then $\bm{\mathrm{r}}_{i-\frac{1}{2},j}=\bm{\mathrm{x}}_{i-\frac{1}{2}}-\bm{\mathrm{x}}_{0}$ is computed and used in the cross product. Mathematically, this is realized as

[TABLE]

in which

[TABLE]

C.3 Change in body momentum

The final term that needs to be discretely approximated is the change in momentum of the immersed body. This is presented as an integral of the body’s velocity over the region $V_{b}(t)\subset\Omega$ in an Eulerian reference frame. Since the body’s position and velocity are described in a Lagrangian reference frame over a region $I_{b}\subset U$, it is generally much easier to evaluate the Lagrangian form of this integral instead. Using the definition of $\delta(\bm{\mathrm{x}})$, it can be shown that the momentum over these two different reference frames are equivalent:

[TABLE]

Letting $\mathcal{G}^{n}$ denote the collection of discrete IB points corresponding to the region $I_{b}$ at time step $n$, the object’s momentum is obtained by

[TABLE]

in which $\bm{\mathrm{U}}_{l,m}^{n}$ denotes the velocity of IB node $(l,m)$ at time step $n$, and $\Delta\bm{\mathrm{s}}_{l,m}$ denotes the discrete volume occupied by the node. The change in momentum required for the evaluation of hydrodynamic forces is then given by

[TABLE]

The change in the body’s angular momentum for the torque calculation is done similarly.