Improving convergence of volume penalised fluid-solid interactions

TL;DR

This paper analyzes and improves the volume-penalty method for fluid-solid interactions, achieving second order accuracy by smoothing the damping and applying Richardson extrapolation, thus enhancing convergence and reducing model error.

Contribution

The authors introduce a smoothing prescription for damping in volume penalization, reducing model error to second order and improving convergence in fluid-solid interaction simulations.

Findings

Displacement length causes dominant model error proportional to Re and damping time.

Smoothing the damping eliminates displacement length, reducing error to O(η).

Combining Richardson extrapolation with smoothing achieves O(η^2) convergence.

Abstract

We analyse and improve the volume-penalty method, a simple and versatile way to model objects in fluid flows. The volume-penalty method is a kind of fictitious-domain method that approximates no-slip boundary conditions with rapid linear damping inside the object. The method can then simulate complex, moving objects in general numerical solvers without specialised algorithms or boundary-conforming grids. Volume penalisation pays for this simplicity by introducing an equation-level error, the $\textit{model error}$, that is related to the damping time $\eta \ll 1$. While the model error has been proven to vanish as the damping time tends to zero, previous work suggests convergence at a slow rate of $\mathcal{O}(\eta^{1/2})$. The stiffness of the damping implies conventional volume penalisation only achieves first order numerical accuracy. We analyse the volume-penalty method using multiple-scales matched-asymptotics with a signed-distance coordinate system valid for arbitrary smooth geometries. We show the dominant model error stems from a displacement length that is proportional to a Reynolds number $\text{Re}$ dependent boundary layer of size $\mathcal{O}(\eta^{1/2}\text{Re}^{-1/2})$. The relative size of the displacement length and damping time leads to multiple error regimes. Our key finding derives a simple smoothing prescription for the damping that eliminates the displacement length and reduces the model error to $\mathcal{O}(\eta)$ in all regimes. This translates to second order numerical accuracy. We validate our findings in several comprehensive benchmark problems and finally combine Richardson extrapolation of the model error with our correction to further improve convergence to $\mathcal{O}(\eta^{2})$.

Tables

Table 1. Table 1: Optimal zero-shift smoothing δ ∗ ​ ( ℓ = 0 ) superscript 𝛿 ℓ 0 \delta^{*}(\ell=0) for different normalised mask profiles Γ ^ ^ Γ \widehat{\Gamma} .

Mask profile	${\hat{\Gamma }}_{\mathrm{erf}}$	${\hat{\Gamma }}_{[\mathrm{erf};1]}$	${\hat{\Gamma }}_{\tanh }$	${\hat{\Gamma }}_{[\tanh ;1]}$
Optimum smoothing	3.11346786	3.80171928	2.64822828	3.54403048

Table 2. Table 2: Table of fixed cylinder simulation parameters.

Parameter	Symbol	Value
Reynolds number	$\mathrm{Re}$	200
Inner domain radius	$R_1$	0.1
Outer domain radius	$R_2$	10
Final time	$t_{end}$	10
Azimuthal modes	$N_{\theta }$	512

Table 3. Table 3: Table of varying cylinder simulation parameters.

Simulation	$\ell$	$\delta$	$\eta$	$\epsilon$
$(0,0)$	0	0	$1\text{×}{10}^{-2}$	$7.07\text{×}{10}^{-3}$
$(0,1)$	0	0	$1\text{×}{10}^{-3}$	$2.23\text{×}{10}^{-3}$
$(0,2)$	0	0	$1\text{×}{10}^{-4}$	$7.07\text{×}{10}^{-4}$
$(1,0)$	$7.07\text{×}{10}^{-2}$	0	$1\text{×}{10}^{-2}$	$7.07\text{×}{10}^{-2}$
$(1,1)$	$7.07\text{×}{10}^{-3}$	0	$1\text{×}{10}^{-3}$	$2.23\text{×}{10}^{-3}$
$(1,2)$	$7.07\text{×}{10}^{-4}$	0	$1\text{×}{10}^{-4}$	$7.07\text{×}{10}^{-4}$
$(2,0)$	0	$2.69\text{×}{10}^{-2}$	$1\text{×}{10}^{-2}$	$7.07\text{×}{10}^{-2}$
$(2,1)$	0	$8.50\text{×}{10}^{-3}$	$1\text{×}{10}^{-3}$	$2.23\text{×}{10}^{-3}$
$(2,2)$	0	$2.69\text{×}{10}^{-3}$	$1\text{×}{10}^{-4}$	$7.07\text{×}{10}^{-4}$

Equations

$${\partial }_{t}u+\frac{1}{\rho }\nabla p-\nu {\nabla }^{2}u+u\cdot \nabla u-f$$

$$u=u_s\text{on}\partial {\Omega }_s,$$

$$F$$

$$\Sigma =-pI+\rho \nu (\nabla u+\nabla u^{\top }).$$

$${\partial }_{t}u+\frac{1}{\rho }\nabla p-\nu {\nabla }^{2}u+u\cdot \nabla u-f$$

$$F$$

$$T$$

$$\mathrm{Re}$$

$$\eta =\mathrm{Re}{\epsilon }^2,$$

$${\partial }_{t}u+u\cdot \nabla u+\nabla p-\frac{1}{\mathrm{Re}}{\nabla }^{2}u-f$$

$$\displaystyle\widehat{\Gamma}^{c}(x)=\widehat{\Gamma}\Big{(}\tfrac{x}{\sqrt{1-{x^{2}}/{c^{2}}}}\Big{)}\quad\text{where}\quad-c<x<c.$$

$${\Gamma }_{\ell ,\delta }(x)$$

$${\hat{\Gamma }}_{\tanh }(x)$$

$${\mathcal{E}}_1=||u-u_0|{|}_{L^1({\Omega }_f)}=\frac{1}{|{\Omega }_f|}{\int }_{{\Omega }_f}|u-u_0|dV.$$

$${\mathcal{E}}_{\infty }=||u-u_0|{|}_{L^{\infty }({\Omega }_f)}=\underset{x\in {\Omega }_f}{\max }|u-u_0|.$$

$$\Delta F$$

$$x=p(s)+\sigma \hat{n}(s).$$

$$t_i$$

$${\nabla }_{\perp }$$

$$\hat{n}$$

$$u=u_{\sigma }\hat{n}+u_{\perp },\text{where}{u}_{\perp }=u_1{\hat{t}}_1+u_2{\hat{t}}_2.$$

$$\nabla$$

$$\overline{K}$$

$$\hat{K}$$

$$\nabla \cdot u$$

$${\nabla }^{\perp }$$

$$u^{\perp }$$

$$\nabla \times u$$

$$-\nabla \times \nabla u$$

$$+{\hat{J}}^{-1}{\partial }_{\sigma }\left(\hat{J}{J}^{-1}({\partial }_{\sigma }(J{u}_{\perp })-{\nabla }_{\perp }{u}_{\sigma })\right)+{\hat{J}}^{-1}{\nabla }^{\perp }\left(\frac{{\nabla }^{\perp }\cdot (J{u}_{\perp })}{|J|}\right).$$

$${\nabla }_{\perp }\hat{n}$$

$${\mathcal{R}}_1^{12}$$

$$\nabla u$$

Full text

1 Introduction

Fluid-solid interactions are a continuing challenge for numerical modelling. The main difficulty results from imposing fluid boundary conditions at the interface of complicated solid shapes. In the most interesting cases, the location, shape, or topology of a solid co-evolves with the dynamics of the fluid.

The most obvious way to handle complex boundary conditions is with “boundary-conforming” grids or numerical elements. While many sophisticated numerical packages take this approach [17, 20, 33, 13, 4], the method is expensive to simulate, and challenging to implement from a software perspective. There is a need for simpler schemes that are compatible with existing codes designed to solve standard fluid-dynamical models in simple domains.

This paper examines the volume-penalty method (VPM), which approximates no-slip boundaries with rapid damping inside the fictitious solid interior. The VPM thereby avoids explicit tracking of the boundary; dynamically updating source terms is comparatively easy. Fluid-stress boundary integrals are also replaced by volume integrals of the damping. This simple method arose from general forcing approaches [7, 25], and is motivated physically by porous-media flows [2, 3]. Its versatility has spawned applications to simulations of the Navier-Stokes equations with arbitrary geometries [35, 34, 29, 52], moving boundaries [42, 36], magnetohydrodynamics [48], flux boundary conditions [31, 46], and even insect flight [18, 19].

Philosophically, boundary conditions exist to model rapid media changes. A “solid” wall is, in reality, an extreme stiffness change in the elastic properties of a continuum. Students in physics often first learn about boundary conditions for Maxwell’s equations with “pill-box” arguments applied to rapidly-changing macroscopic media. The VPM simply backs off from these notions and allows some level of resolved rapid change in medium properties.

There is a long history of this approach, disguised in many forms. The Rankine-Hugoniot jump conditions emerge in the limit of vanishing diffusivity in Burgers equation [16, 28]. Gibbs’ analysis of capillary forces [22] and Stefan’s melting phase transitions [49] supposed material discontinuities. Further studies found small-but-finite length scales standing in for infinitesimal jumps (e.g., Van der Waals’ explanation of surface tension [51], and Cahn and Hilliard’s model of phase separation [12]).

Many methods exploit smoothness to design simple algorithms for simulating PDEs on complicated domains. Phase-field models are a prime example of this smooth approach [27, 38, 5]. Diffuse-interface methods extend the phase-field approach to simulate two-phase liquids [30, 1]. The diffuse-domain method also originates from phase-field models and approximates Neumann, Robin, and Dirichlet boundary conditions [40]. Fictitious-domain methods are similar to diffuse-domain methods and solve for Lagrange multipliers within the object to enforce the desired boundary conditions [24, 23, 55, 54]. The immersed-boundary method is another classical method for simulating elastic objects in fluids [43, 44] (though rigid objects can cause numerical stiffness [42]). Of course, the length scales of realistic medium transitions are often too small to resolve computationally. In these situations, smooth models are a numerical appliance which trade the accuracy of boundary conditions for greater computational simplicity.

The key challenge of the VPM, and other smoothed models, is to minimise the error compared with formal boundary conditions. Much progress has been made in improving the error of phase-field models [10, 11, 32, 15] and diffuse-domain methods [21, 39, 8, 45, 53]. Progress has been slower for the VPM.

Previous work has shown that the error of the VPM vanishes in the limit that the damping timescale $\eta$ tends to zero [2]. Later work showed this error converges at the slow rate $\mathcal{O}(\eta^{1/2})$ [14]. The stiffness of the damping implies this error leads to merely first-order numerical convergence. These works also used discontinuous damping, which restricts numerical convergence, causing Gibbs phenomena in spectral codes [34].

In this paper, we improve upon previous understanding of the VPM and derive simple corrections that reduce the overall error. Our analysis extends the higher-order analysis of phase-field models [32, 15] and diffuse-domain methods [45, 55]. We begin in section 2 by describing the VPM, and the various metrics used to measure the error. The central penalty quantity is the non-dimensional damping length scale, $\varepsilon$, which depends on the Reynolds number Re and non-dimensional damping time scale, $\eta$. Section 3 then develops a general and straightforward multiple-scales matched-asymptotics procedure in powers of $\varepsilon$. We use this analysis to quantify the error for arbitrary smooth objects in three dimensions. In section 3.1, we introduce the signed-distance coordinate system. This analytical tool helps streamline complicated calculations and unambiguously accounts for perplexities such as interface curvature. We analyse the VPM up to second-order and derive simple conditions to eliminate the leading-order error. Section 4.1 develops efficient methods to calculate these optimisations in general cases. The remainder of the paper verifies the improvement to performance in several test problems. In section 5, we test our prescriptions in time-dependent accelerating flow past a rotating 2D cylinder. We outline the efficient spectral code Dedalus, use it to generate numerical reference and penalised solutions, and validate the improved accuracy of our corrections. We finally combine our optimisations with Richardson extrapolation of the model error to further reduce model error to $\mathcal{O}(\eta^{2})$.

2 The volume-penalty method

2.1 Objects in incompressible hydrodynamics

Fluid-solid interactions in incompressible viscous fluids are normally modelled by partitioning the domain $\Omega$ into solid $\Omega_{s}$ and fluid $\Omega_{f}$ components (see fig. 1 $(a)$). On $\Omega_{f}$ we solve the Navier-Stokes equations,

[TABLE]

for the fluid velocity ${u}$ and the pressure $p$, given a kinematic viscosity $\nu$, constant density $\rho$, and external body force ${f}$. The solid then acts through no-slip boundary conditions

[TABLE]

where ${u}_{s}$ is the solid velocity at the boundary. The fluid in turn influences the solid through drag ${F}$ and torque ${T}$ surface stress integrals

[TABLE]

where $\widehat{n}$ is the outward unit normal vector on the solid boundary $\partial\Omega_{s}$, and $\Sigma$ is the stress

[TABLE]

Unfortunately this formulation can be difficult to discretise in elaborate geometries. This motivates the volume-penalty method, an alternate model of fluid-solid interactions which allows simpler discretisation, at the cost of some error in the mathematical solution.

2.2 The volume-penalty approximation

The volume-penalty method models objects in fluids by replacing no-slip boundary conditions with an added linear damping in the momentum equation,

[TABLE]

where the mask function $\Gamma:\Omega\to[0,1]$ represents the location of the solid (which may evolve in time). The penalty term has a physical interpretation. It corresponds to the Darcy drag of the Brinkman equations for fluid flow in a porous medium [3]. The small damping timescale $\tau$ forces the fluid velocity to tend to the solid velocity within the object. Traditionally, the penalty mask $\Gamma$ is chosen to be a discontinuous indicator function over the solid region [2, 3, 37, 48, 52, 47, 18, 19]. However, we emphasise that this is not the only choice. We will later show that it is in fact suboptimal.

The volume-penalty method also approximates the force $F$ and torque $T$ surface integrals (eq. 3, [7]) with simpler-to-discretise volume-penalty integrals [2]

[TABLE]

The second terms account for forcing within the object $\Omega_{s}$ and are analogous to the unsteady correction from [18, 50]. The final corrections account for possible no-slip boundaries $\partial\Omega_{in}$ within the volume penalised object (as in fig. 1). Both corrections are necessary in section 5.

The volume-penalty method is useful because it removes the object boundary from the mathematical formulation. A single set of equations applies throughout the entire domain, the solid is represented by a mask function $\Gamma$, and the force and torque are calculated with volume integrals. This makes it simple to implement numerically, as we do not need the numerical grid to conform to the solid, nor do we need to interpolate quantities at the boundary.

However, the volume-penalty model is different to the standard no-slip model. The mathematical solution of each problem will therefore differ. We define this equation-level difference as the model error.

Definition 2.1** (Model error).**

The difference between the mathematical solution of the true no-slip boundary formulation, and the volume-penalty approximation. This excludes numerical differences due to discretisation.

We emphasise that the model error exists separate to the discretisation error, which is specific to a particular numerical implementation. The model error is implementation-independent, and therefore more fundamental. By finding mathematical ways to minimise the model error, we derive improvements that apply to all possible implementations of the volume-penalty method.

But what is the model error? And what should it depend on? We answer these questions in the remainder of this section.

2.3 Non-dimensional penalty parameters

The most important determinants of the model error are two non-dimensional penalty parameters. We initiate our non-dimensionalisation by defining the dimensional damping time $\tau$, the characteristic length scale of the domain $L$, and the characteristic velocity scale of the fluid $U$. These induce three key non-dimensional parameters, the Reynolds number $\operatorname{Re}$, the damping time scale $\eta$, and the damping length scale $\varepsilon$

[TABLE]

The Reynolds number $\operatorname{Re}$ characterises the ratio of viscous and inertial time scales, the (non-dimensional) damping time scale $\eta$ measures the separation of damping and inertial time scales, and the damping length scale $\varepsilon$ comes from balancing the viscous and damping terms.

The damping time $\eta$ and length $\varepsilon$ are closely related via the Reynolds number $\operatorname{Re}$

[TABLE]

and are equal when $\eta=\operatorname{Re}^{-1}=\varepsilon$. The length and time scale each induce different kinds of model errors. Either side of this equality one scale dominates, leading to distinct error regimes.

We then write the penalised equations as

[TABLE]

where $u$ and $p$ are penalised solutions. Solutions to the ‘true’ non-dimensional Navier-Stokes equations with no-slip boundary conditions are denoted with a zero subscript (e.g. $u_{0}$, $p_{0}$), in keeping with the convergence of the penalised solution $u$ to the true solution $u_{0}$ as $\varepsilon\to 0$.

Remark 2.1** (Error vs computational effort).**

To numerically resolve the damping length $\varepsilon$ and time scale $\eta=\operatorname{Re}\varepsilon^{2}$, the spatial degrees of freedom (Fourier modes, grid points etc.) must scale as $\varepsilon^{-D}$ in dimension $D$, and total time steps must scale as $\varepsilon^{-2}$, for total effort $\varepsilon^{-(D+2)}$. Given an error scaling of $\varepsilon^{\alpha}$, the increased effort to halve the error is proportional to $2^{(D+2)/\alpha}$. This implies a 16-fold increase in cost in 2D and 32-fold increase in 3D for the existing scaling $\alpha=1$. It is thus essential to improve the error scaling $\alpha$.

2.4 Mask functions

The final quantity that influences the model error is the choice of mask function $\Gamma$. We now define normalised mask functions $\widehat{\Gamma}$, and the general masks $\Gamma$ that we consider in this paper.

Definition 2.2** (Normalised mask function $\widehat{\Gamma}$).**

A function on $\mathbb{R}$ which satisfies

Boundedness: $\widehat{\Gamma}:\mathbb{R}\to[0,1]$. 2. 2.

Limiting behaviour: $\lim_{x\to-\infty}\widehat{\Gamma}(x)=1$ and $\lim_{x\to\infty}\widehat{\Gamma}(x)=0$. 3. 3.

Monotonicity: $x_{2}>x_{1}\implies\widehat{\Gamma}(x_{2})<\widehat{\Gamma}(x_{1})$. 4. 4.

Symmetry: $\widehat{\Gamma}(x)+\widehat{\Gamma}(-x)=1$. 5. 5.

Normalisation: \frac{d\widehat{\Gamma}}{dx}\big{|}_{x=0}=-1.

We concentrate on smooth normalised mask functions and more restrictive compact normalised mask functions. Any smooth normalised mask function $\widehat{\Gamma}$ can be compactified to the interval $[-c,c]$ using the following transformation [6],

[TABLE]

General mask functions $\Gamma_{\ell,\delta}$ can be generated from a normalised mask function $\widehat{\Gamma}$ through scaling $\delta$ or shifting $\ell$,

[TABLE]

Discontinuous masks are the limiting case of smooth functions as the thickness approaches zero. The standard approach takes the limit $\delta\to 0,\ell=0$ first, then considers convergence as $\varepsilon\to 0$. This is not optimal. Choosing a particular $\ell$ or $\delta\propto\varepsilon$ cancels the leading order error of the volume-penalty method.

For concreteness, we define the hyperbolic tangent mask $\widehat{\Gamma}_{\tanh}$ and error function mask $\widehat{\Gamma}_{\operatorname{erf}}$,

[TABLE]

and compactified versions thereof, $\widehat{K}_{[\tanh;c]},\widehat{K}_{[\operatorname{erf};c]}$, on the interval $[-c,c]$.

These definitions can be extended to higher dimensions by defining a locally conformal coordinate system at the mask boundary and applying the one dimensional mask functions in the wall-normal direction. Before developing this coordinate system in section 3.1, we now give several metrics of the model error.

2.5 Model error metrics and the displacement length

There are many possible ways to quantify the model error of the volume-penalty method. We give two metrics that measure the accuracy in an average and worst-case sense, as well as the physically meaningful error in the drag calculation.

Definition 2.3** (Global error $\mathcal{E}_{1}$).**

The $L^{1}$ norm of the difference of the penalised solution $u$ from the reference solution $u_{0}$ in the fluid domain,

[TABLE]

Definition 2.4** (Local error $\mathcal{E}_{\infty}$).**

The $L^{\infty}$ norm of the difference of the penalised solution $u$ from the reference solution $u_{0}$ in the fluid domain,

[TABLE]

Definition 2.5** (Drag error $\Delta F$, Torque error $\Delta T$).**

The difference between the volume integral for the penalised drag $F$ (torque $T$) and the surface stress integral for the reference drag $F_{0}$ (reference torque $T_{0}$) on a true no slip boundary.

[TABLE]

We also justify the fluid error norms by noting that together they quantify error localisation. If $\left|u-u_{0}\right|\approx\varepsilon$ in a region of size $d$, and much smaller elsewhere, then $\mathcal{E}_{1}\approx\varepsilon d$ (for a unit volume domain), while $\mathcal{E}_{\infty}\approx\varepsilon$, so the ratio $\mathcal{E}_{1}/\mathcal{E}_{\infty}\approx d$ quantifies the degree of error localisation.

The key insight of this paper is that conventional masks do not optimise the global error $\mathcal{E}_{1}$ or the drag error $\Delta F$ with respect to the penalty parameter $\varepsilon$. This is because in general the far field behaviour of the penalised velocity $u$ is shifted by an amount proportional to the damping length scale $\varepsilon$. We call this difference the displacement length $\ell^{*}$,

Definition 2.6** (Displacement length $\ell^{*}$).**

(i). The difference in size between the desired no-slip boundary $\partial\Omega_{s}$, and a hypothetical no-slip boundary $\partial\Omega_{s}^{*}$ the mask function $\Gamma$ “most closely” approximates. (ii). Alternatively, the difference in size between the desired no-slip boundary $\partial\Omega_{s}$, and the optimal mask function $\Gamma^{*}$ that most closely approximates it.

This definition of the displacement length can refer to the size of the mask or the desired no-slip boundary. We use these meanings interchangeably as they relate to the same phenomenon: the ideal mask is in general not the same shape or size as the solid we wish to approximate. The concrete advance of this paper is a simple correction to offset the displacement length error by refining the edge of the mask. This localises the error to the boundary layer, reducing the average fluid error $\mathcal{E}_{1}$ and the drag and torque errors $\Delta F$ and $\Delta T$ to $\mathcal{O}(\varepsilon^{2})$. Having defined the appropriate concepts by which to measure the model error, we now analyse it using a general multiple-scales matched-asymptotics framework.

3 Multiple-scales matched-asymptotics analysis of the model error

Understanding volume penalisation requires analysis of the model error. This analysis must be performed at the equation level, and should be a straightforward procedure that works in general geometries and flow configurations. We provide such a framework using multiple-scales matched-asymptotics valid for general tensor-valued equations in arbitrary smooth geometries. This procedure breaks into several steps:

Partition the domain into solid $\Omega^{-}$, fluid $\Omega^{+}$, and interfacial regions $\Delta\Omega$ (fig. 2 $(a)$). 2. 2.

Use the signed distance $\sigma$ as a coordinate in the interfacial region (section 3.1, fig. 2 $(b)$). 3. 3.

Rescale the normal coordinate and operators by $\varepsilon$ in the interfacial region (section 3.2). 4. 4.

Expand the variables in a formal asymptotic power series in $\varepsilon$ each region (section 3.3). 5. 5.

Connect regions with and asymptotic matching conditions (section 3.4). 6. 6.

Iterate to solve the zeroth, first, and second order problems (sections 3.5, 3.6 and 3.7).

The philosophy of the approach is to determine the asymptotic evolution of the volume-penalty method accurate to second order in $\varepsilon$. We show that it is possible to analytically determine optimal damping and mask parameters which eliminate the order $\varepsilon$ model error of the volume-penalty method. We perform the analysis in a body-centred frame incorporating possible acceleration terms into $f$.

The first step partitions the domain into the fluid $\Omega^{+}$, solid $\Omega^{-}$, and interfacial regions $\Delta\Omega$, distinguished by the behaviour of the mask function $\Gamma$. Understanding the interfacial region calls for an appropriate choice of coordinates.

3.1 Signed-distance coordinates

To understand the volume-penalty method we must analyse the boundary region. An ideal coordinate system for boundary layers uses the signed distance $\sigma$ of each point to the closest point on the surface. This choice conforms to the object, aligning the basis vectors with the boundary, and is orthogonal, simplifying the vector calculus.

The signed distance $\sigma$ is defined as the distance of a point $x$ to its nearest point on the interface $p$. It is not difficult to show the point $x$ must also lie in the direction of the unit normal vector $\widehat{n}$ from point $p$, which we label with surface coordinates $s$ (fig. 2 $(b)$),

[TABLE]

The chain rule allows us to rewrite the normal and surface partial derivatives in terms of the gradient operator, $\partial/\partial\sigma=\widehat{n}\cdot\nabla$, $\partial/\partial s_{i}=(\partial p/{\partial}s_{i})\cdot(I+\sigma\left.\nabla\widehat{n}\right|_{\sigma=0})\cdot\nabla$. While the normal derivative is straightforward, the tangential derivatives require an understanding of surface differential geometry. Only then can we invert to find the gradient in the new coordinates.

The first geometric quantity induced by our choice of surface coordinates is a tangent vector basis $t_{i}$. Presuming orthogonal surface coordinates, we also derive the dual vector basis $\nabla s_{i}$ and a unique orthonormal tangent basis $\widehat{t}_{i}$,

[TABLE]

This gives us the surface gradient $\nabla_{\bot}$. Integrating by parts using the surface area measure $dA=|t_{1}||t_{2}|$ gives us the surface divergence $\nabla_{\bot}\cdot$ of a tangent vector field $u_{\bot}=u_{1}\widehat{t}_{1}+u_{2}\widehat{t}_{2}$,

[TABLE]

It is not difficult to show the normal $\widehat{n}$ is everywhere equal to the signed distance gradient. This implies the gradient of the normal is symmetric, and hence diagonalisable. The eigenvectors are the principal directions of curvature (and must align with orthogonal surface coordinates), and the eigenvalues are the principal curvatures $\kappa_{i}$,

[TABLE]

With an orthonormal frame defined throughout the interfacial region, we decompose the fluid velocity into normal and tangential directions,

[TABLE]

We can now invert all partial derivative operators and re-express the gradient

[TABLE]

where we have defined the scale factor $J$, which possesses several useful identities relating to the mean curvature $\overline{K}$, the Gaussian curvature $\left|K\right|$, and volume measure $dV$,

[TABLE]

Using the adjugate tensors $\widehat{J}$ and $\widehat{K}$, we determine the divergence and scalar Laplacian,

[TABLE]

By defining the cross product with the unit normal, which admits simple identities,

[TABLE]

we can take the curl, and hence find the divergence-free vector Laplacian,

[TABLE]

We next find the gradient of the basis vectors, and define Ricci rotation coefficients,

[TABLE]

which finally gives us the gradient and advective derivatives of general vector fields

[TABLE]

3.2 Rescaling interfacial coordinates

The interfacial region is defined as the size $\varepsilon$ region in which the fluid transitions to the solid. Analysing this region requires rescaling the normal coordinate by $\varepsilon$. This in turn projects each differential operator into a formal asymptotic series as follows,

[TABLE]

[TABLE]

3.3 Variable expansions with formal power series

After splitting the domain into the fluid $\Omega^{+}$, solid $\Omega^{-}$, and interfacial $\Delta\Omega$ regions, each variable $q$ in each region (fluid $q^{+}$, solid $q^{-}$, and interfacial $q$) is expressed as a power series in $\varepsilon$,

[TABLE]

We perform this expansion for the fluid velocity $u=u_{\sigma}\widehat{n}+u_{\bot}$ and pressure $p$. We then substitute these variable expansions into the infinite hierarchy of equations generated by section 3.2 to derive a system of equations at each order of $\varepsilon$. Collecting terms and solving the model equations at each order requires a matching procedure between adjacent regions, which we describe in section 3.4.

3.4 Asymptotic matching

To ensure agreement between different regions we specify asymptotic matching boundary conditions. This subtle notion requires asymptotic agreement in intermediate zones $\xi\sim\varepsilon^{-1/2}$ in the limit that $\varepsilon\to 0$. We can then let $\xi$ approach infinity for the inner problem without encountering coordinate singularities (provided $\varepsilon\ll\min_{i}|\kappa_{i}^{-1}|$), and let $\sigma$ approach zero for the outer problem without entering the interfacial region. That is, for any variable $q$, we require

[TABLE]

Each variable is expressed as an asymptotic series. Each series term in the outer variables can be further expanded as a Taylor series about the interface, involving normal derivatives $\partial_{\sigma}$ and the rescaled signed distance $\varepsilon\xi=\sigma$. Combining these expansions and solving order by order in $\varepsilon$, the matching conditions at each order simplify to

[TABLE]

Written explicitly, we expect the limiting behaviour of the inner expansions to be asymptotically equivalent to interfacial boundary conditions on the outer variables as $\xi\to\pm\infty$,

[TABLE]

Equipped with our general multiple scales matched asymptotic procedure for arbitrary geometries, we can now apply it to the volume-penalty method and derive its asymptotic behaviour in a straightforward step-by-step manner.

3.5 Zeroth order: The tangential boundary layer

Fluid problem — At leading order in $\varepsilon$ we reproduce the fluid equations eq. 1

[TABLE]

Solid problem — The leading order solid problem is

[TABLE]

Boundary problem — The leading order boundary problem illustrates the structure of the subsequent orders. The divergence condition reduces to purely the normal direction, and implies the normal velocity is constant. Matching to the solid problem requires this to be zero, giving a no-flux boundary condition for the fluid velocity,

[TABLE]

Tangential velocity boundary layer — We now analyse the tangential velocities with the leading order momentum equation, which takes the form

[TABLE]

We emphasise this leading order operator, as it is key to understanding the volume-penalty method. Noting the asymptotic behaviour of $\Gamma$, we can infer the asymptotic behaviour of the kernel. The first solution $\mathcal{V}$ grows unphysically in the solid, and is constant in the fluid

[TABLE]

The physical solution $\mathcal{U}$ decays exponentially into the solid, and grows linearly in the fluid

[TABLE]

The shift $\ell^{*}$ is the displacement length of the mask. Example solutions for different $\Gamma$ are shown in fig. 4. The solid matching condition tells us that $u_{\bot 0}$ is some multiple of the $\mathcal{U}$ solution,

[TABLE]

The fluid matching condition instead requires zero gradient in the fluid. This is only satisfied by the zero solution, implying no-slip boundary conditions at $\partial\Omega_{s}$ for the fluid velocities,

[TABLE]

The leading order problem is therefore equivalent to the incompressible Navier-Stokes equations with no-slip boundaries at $\sigma=0$, justifying convergence of the volume-penalty method. We continue with the next order problem, the leading order error, and show how it can be controlled through judicious choice of $\Gamma$.

3.6 First order: Displacement length error

Fluid problem — In the fluid we now have

[TABLE]

Solid problem — The solid then gives

[TABLE]

Boundary problem — The continuity equation for the boundary problem implies the first order normal velocity is again constant, which is zero by matching with the solid boundary condition

[TABLE]

Hence the normal velocity is at most $\mathcal{O}(\varepsilon^{2})$ around the boundary. The normal component of the momentum equation implies the leading order behaviour of the pressure is constant

[TABLE]

Error for general masks — The tangential ($\bot$) momentum equation implies

[TABLE]

Matching to the solid problem implies $u_{\bot 1}$ is some multiple of $\mathcal{U}$

[TABLE]

and matching the fluid problem gives us the boundary condition for $u_{\bot 1}^{+}$,

[TABLE]

This shift $\ell^{*}$ is the displacement length of the mask, and depends only on $\Gamma$. Most masks give displacement length $\ell^{*}$ order unity (see section 4.1 for examples), implying the boundary condition $u^{+}_{\bot 1}(0)=\mathcal{O}(\varepsilon)$. This boundary perturbation of the fluid problem causes first order error $\mathcal{O}(\varepsilon)$ throughout the fluid.

**Error for optimal masks ** — However, if we choose a mask with zero displacement length ($\ell^{*}=0$), we obtain homogenous boundary conditions in the fluid,

[TABLE]

If the errors $u^{+}_{1}$ and $p^{+}_{1}$ are initially zero, the homogeneous linear momentum equation implies they will remain so for all time. In reality, any small perturbation will grow over time (according to the Lyapunov time scale of the flow), but this is inherent to the chaotic nature of most fluid-flows, and not a limitation of this specific method. The statistical behaviour of the fluid will be accurate to $\mathcal{O}(\varepsilon^{2})$ with this correction. We show how to choose a mask with no displacement error in section 4.1, but first analyse the second order problem, and in doing so find the emergence of two $\operatorname{Re}$ dependent asymptotic regimes.

3.7 Second order: Different volume-penalty regimes

Fluid problem — We have an added inhomogeneity from the nonlinear first order convection term, which is zero for optimised masks,

[TABLE]

Solid problem — We now have non-zero flow due to pressure gradients and forcing in the solid,

[TABLE]

Boundary problem — The incompressibility equation simplifies since $u_{\sigma 1}=0$, giving $u_{\sigma 2}$

[TABLE]

We can simplify noting that at the boundary $\sigma=0^{\pm}$,

[TABLE]

where the first equality comes from matching to the solid solution, and the second equality comes from evaluating the normal component of the fluid momentum equation at $\sigma=0$. This behaves as a constant plus an exponential in the solid, and quadratically within the fluid,

[TABLE]

The quadratic behaviour is consistent with the leading order problem, and the linear behaviour matches onto zero for an optimised mask ($A_{\sigma 2}\propto\ell$). However, the constant behaviour into the fluid depends on the jump in the pressure derivative (and external forcing) across the boundary, which depends on the details of the flow, and cannot be eliminated using a passive mask, so in general $B_{\sigma 2}$ is non-zero. Hence, the boundary condition for the second order normal velocity in the fluid scales with the Reynolds number $\operatorname{Re}$ (for a given jump in the pressure).

We can also take the divergence of the momentum equation in the boundary region to obtain

[TABLE]

Integration shows a discontinuity in the normal pressure gradient across the mask boundary.

We can then determine the tangential velocity using the tangential momentum equation,

[TABLE]

This is solved using variation of parameters,

[TABLE]

where $\mathcal{R}$ is the right hand side of eq. 41 and $c$ is a constant of integration.

Reynolds dependent error regimes — The second order error in the fluid scales with the Reynolds number $\operatorname{Re}$, and so is proportional to $\eta$, while the first order error scales with $\varepsilon$. For sufficiently large $\operatorname{Re}>\varepsilon^{-1}$, the asymptotic separation is invalid, and the “second order” time scale error $\mathcal{O}(\operatorname{Re}\,\varepsilon^{2})=\mathcal{O}(\eta)$ dominates the “first order” length scale error $\mathcal{O}(\varepsilon)$. This leads to several asymptotic damping regimes. We summarise these regimes in terms of $\log\eta$ and $\log\operatorname{Re}$ in fig. 3.

If $\eta>1$ (the left half plane), the fluid is damped slowly and the velocity in the solid is large. Likewise if $\varepsilon>1$ (the lower left diagonal half-plane), damping cannot overcome the fluid viscosity. It is possible to have either $\varepsilon$ or $\eta$ be small and the other large, but both must be small to approximate a solid.

Ignoring these large error regimes, we can have intermediate damping (light gray) when the Reynolds number is large enough to promote the second order system to dominance ($1>\eta>\varepsilon>\operatorname{Re}^{-1}$). This regime gives the appearance of $\mathcal{O}(\varepsilon^{2})$ convergence even for unoptimised masks, but this is only temporary, as seen in [3].

Once the damping is sufficiently large ($1>\operatorname{Re}^{-1}>\varepsilon>\eta$) the damping length scale $\varepsilon$ always dominates the time scale $\eta$. This is the strong damping regime (white), for which it is necessary to apply mask corrections to promote the convergence to $\mathcal{O}(\varepsilon^{2})$.

Reynolds dependent cost regimes — We can analyse the computational cost regimes for large $\operatorname{Re}$ by comparing the damping length scale $\varepsilon$ and time scale $\eta$ with the intrinsic Kolmogorov length scale $\lambda$ and CFL time scale $t$ in three dimensional turbulence,

[TABLE]

The grid resolution $\Delta x$ and time step size $\Delta t$ must resolve the smallest scale,

[TABLE]

The computational cost $C$ is proportional to the degrees of freedom $\Delta x^{-3}\Delta t^{-1}$. Scaling the damping in terms of the Reynolds number $\eta=\operatorname{Re}^{-\alpha}$, we find several cost regimes for $\operatorname{Re}\gg 1$,

[TABLE]

Volume penalised boundaries can thus be added at no extra cost with an error of $\operatorname{Re}^{-1/2}$. Beyond this regime $\alpha>1/2$, the computational cost scales more aggressively in $\eta$. This has given the impression that low $\operatorname{Re}$ simulations are difficult using volume penalisation. However, smaller $\operatorname{Re}$ makes it easier to achieve strong damping where mask corrections can now allow $\mathcal{O}(\varepsilon^{2})$ convergence. We therefore extend the practicability of the volume-penalty method to low $\operatorname{Re}$. We finally note that this analysis is premised on being in the strong damping regime, where the second order error is smaller than the first order error.

3.8 Force and torque calculations

We now relate the integral of the penalty term to the force and torque calculations. For a general mask, this integral can be rewritten using the penalised Navier-Stokes equations

[TABLE]

We use a body-centered coordinate frame and account for the acceleration terms in the body forcing $f$. For a mask with support on $\Omega\setminus\Omega^{+}$, which is $\mathcal{O}(\varepsilon)$ larger than $\Omega_{s}$, we can restrict the integral to this region, with exponentially small error for non-compact masks like $\Gamma_{\tanh}$.

[TABLE]

The divergence theorem relates this to the stress along $\partial\Omega^{+}$, plus corrections

[TABLE]

We can compare the penalised stress $\Sigma$ with the true stress $\Sigma_{0}$

[TABLE]

We can then relate the true stress on the mask support boundary $\partial\Omega^{+}$ to the true stress on the true boundary $\partial\Omega_{s}$ by Stokes’ theorem with the unpenalised equation for $u_{0},p_{0}$

[TABLE]

Both the stress on interior no-slip boundaries (as in section 5) and the $f$ volume integral are at order $\varepsilon^{0}$, so must be subtracted off (note that $f$ is integrated over exactly $\Omega_{s}$). The leading stress error $\Sigma-\Sigma_{0}$ is $\mathcal{O}(\varepsilon)$ at $\Omega^{+}$ for general masks, as both $p^{+}_{1}$ and $u^{+}_{1}$ are non-zero in the fluid, implying $\mathcal{O}(\varepsilon)$ convergence of the force calculation.

However, an optimal mask ensures that $p^{+}_{1}$ and $u^{+}_{1}$, and hence $\Sigma-\Sigma_{0}$ are accurate to $\mathcal{O}(\varepsilon^{2})$ at the mask boundary $\partial\Omega^{+}$, so this term is now $\mathcal{O}(\varepsilon^{2})$ accurate. The remaining acceleration terms are $\mathcal{O}(\varepsilon^{2})$ in general as both $u_{0}$ and $\varepsilon u_{1}$ are $\mathcal{O}(\varepsilon)$ in a region $\Omega^{+}-\Omega_{s}$ with support $\mathcal{O}(\varepsilon)$. The $\mathcal{O}(\varepsilon)$ error $u_{1}$ decays exponentially fast into the solid, so also has support and magnitude of $\mathcal{O}(\varepsilon)$, integrating to $\mathcal{O}(\varepsilon^{2})$. The $\mathcal{O}(\varepsilon^{2})$ error $u^{-}_{2}$ and $p^{-}_{2}$ have support throughout the solid region, but again integrate to only $\mathcal{O}(\varepsilon^{2})$. Hence, by choosing an optimal mask, the dominant stress error can be eliminated and the force calculation accuracy promoted to $O(\varepsilon^{2})$.

The torque calculation proceeds similarly. We recover the correct stress at $\partial\Omega^{+}$ to $\mathcal{O}(\varepsilon^{2})$, and we need to show this is only $\mathcal{O}(\varepsilon^{2})$ different from the true stress at $\partial\Omega_{s}$. This follows from the divergence theorem applied to $r\times(\nabla\cdot\Sigma_{0})$, which relies on the fact that the stress tensor is symmetric, and that the gradient of the coordinate vector is the identity tensor. Following the same procedure as for the force calculation, one finds $\mathcal{O}(\varepsilon)$ accuracy for the torque for general masks, and $\mathcal{O}(\varepsilon^{2})$ accuracy for optimal masks.

4 Testing optimal masks

In the previous section, we derived the necessary conditions for optimal masks. We now investigate these conditions, and test the performance of our corrections in a series of increasingly general test problems.

4.1 Calculating optimal masks

We first examine the construction of optimal masks for the inner problem of section 3 eq. 22,

[TABLE]

We plot several solutions $\mathcal{U}(\xi)$ to various mask function choices $\Gamma$ in fig. 4.

4.1.1 Conventional volume penalisation

Conventional masks in eq. 22 lead to non-zero displacement length,

[TABLE]

From section 3 the non-zero displacement length $\ell^{*}=-1$ induces an $\mathcal{O}(\varepsilon)$ error in the fluid solution, agreeing with previous estimates of the accuracy of the volume-penalty method. We now examine mask corrections to eliminate this displacement and improve convergence.

4.1.2 Shifting the mask

The simplest correction is to shift the mask distance $\ell=1$ into the fluid (fig. 4 ($a$) vs fig. 4 ($b$)).

[TABLE]

This simple refinement eliminates the displacement error, and section 3 shows this improves the global accuracy $\mathcal{E}_{1}$ and drag and torque errors $\Delta F$ and $\Delta T$ to $\mathcal{O}(\varepsilon^{2})$.

4.1.3 Smoothing the mask

We now demonstrate a key finding of the paper: there exist smooth masks that eliminate the displacement length error. Late in the review stage we became aware of work by Poulsen and Voorhes which independently performed a similar analysis [45]. We consider a $\tanh$ type mask,

[TABLE]

We define transformed coordinates and variables and a rescaled smoothness $n$

[TABLE]

for which the problem becomes

[TABLE]

The solution which satisfies $U(0)=0$ can be written as a Frobenius series

[TABLE]

where here $\Gamma(n)=(n-1)!$ represents the actual gamma function, $\mathrm{B}$ is the beta function, and ${}_{2}F_{1}$ is the hypergeometric function. Considering the behaviour of the function as $z\to 1$ from below, and applying the transformation, we can show that as $\xi\to\infty$

[TABLE]

where $\gamma$ is the Euler–Mascheroni constant and $\psi(n)$ is the digamma function $\Gamma^{\prime}(n)/\Gamma(n)$. We can calculate the value of $n$ (and hence $\delta$) for which the displacement length $\ell(n)$ is zero

[TABLE]

This is the optimal smoothing for an unshifted tanh mask which cancels the displacement error.

We now extend this work to a simple calculation for any mask profile.

4.1.4 General corrections via Riccati transform

We now calculate optimal smoothing and shifting for general mask functions. In current form this requires solving second order boundary value problems. We greatly simplify by using the shifted and scaled coordinate $\zeta$. This transforms the velocity $\mathcal{U}(\xi)$ to the velocity $u(\zeta)$, and a general mask $\Gamma_{\ell,\delta}(\xi)$ to a normalised mask $\widehat{\Gamma}(\zeta)$. We then apply a Riccati transform $R$ to the rescaled eq. 22 to derive an initial value problem for $R(\zeta)$,

[TABLE]

For a mask equal to one for $\zeta<-c$ we have an exponential solution $\mathcal{U}(\xi)\propto e^{\xi}$. Differentiating gives an initial condition for the Riccati equation at $\zeta=c$,

[TABLE]

Solving forward to $R(\zeta=c)$, we then enforce the optimal linear profile $u(\zeta)=\zeta$ beyond $\zeta>c$ to determine the ideal shift $\ell^{*}$ in terms of smoothing $\delta$,

[TABLE]

This transformed system is much easier to solve than the original second order boundary value problem. We solve it using simple Runge-Kutta integrators of the scipy python library. We can then extend to noncompact smooth functions by choosing $c$ sufficiently large.

To illustrate the optimal parameters we plot the (scaled) shift $\ell^{*}$ as a function of smoothing $\delta$ for four normalised masks in fig. 5. The discontinuous limit $\delta\to 0$ reproduces the optimum discontinuous shift $\ell^{*}=1$, and as the smoothing increases, the optimum shift decreases.

An important case is the ‘zero-shift’ optimum smoothing. These are the $\delta^{*}$ intercepts in fig. 5, given in table 1. This choice is numerically motivated by balancing requirements of resolving the mask and the velocity field. Moving away from the optimum curve, the mask is either too ‘large’ ($\ell>\ell^{*}(\delta)$), or too small ($\ell<\ell^{*}(\delta)$).

4.2 Couette flow: Optimal parameters to minimise total error

We have outlined the two sources of error of the volume-penalty method, mathematical model error, and numerical discretisation error. The question faced by practitioners is, given some fixed computational resources, how can one minimise the total error? Damping influences the model error and discretisation error in contrasting ways. Weak damping means the body does not approximate a solid, worsening model error. Strong damping induces the small length scales that worsen the discretisation error if not resolved. There must exist some optimum between these extremes which minimises the sum of the model and discretisation error.

We examine the numerical performance of our corrections in a specific test problem: a Chebyshev discretisation of Couette flow. The mathematical problem is a rescaling of the inner problem which solves for the tangential velocity $v$ near a volume-penalised boundary with a prescribed unit velocity gradient into the fluid

[TABLE]

We solve this problem numerically using Chebyshev polynomials in Dedalus [9] for five choices of numerical resolution, $n=128,192,256,384,$ and $512$ Chebyshev modes. For each resolution $n$ we compare three mask profile strategies: a standard discontinuous mask, an optimal shifted discontinuous mask $\ell=\varepsilon$, and an optimal unshifted smoothed mask $\delta=\delta^{*}\varepsilon$. We choose $\varepsilon=4/n$ to balance the model and discretisation errors. We then measure the total error as the $L^{1}$ and $L^{\infty}$ norm of the difference between the numerical solution and the ideal solution $v_{0}(x)=x$ in the fluid domain, on a resolved grid of 16 times higher resolution than the grid used to calculate the solution.

We plot the results in fig. 6, which shows that smoothed masks achieve lower total error (in $L^{1}$ norm) than any discontinuous mask, for a given numerical resolution. Shifted masks achieve slightly better $L^{\infty}$ error, but at the cost of large Gibbs oscillations that reduce numerical convergence to first order in the fluid domain. We caution that the exact errors are dependent on different numerical discretisations. However, the qualitative findings are general: optimal smoothed masks improve the volume penalty method to second order numerical convergence.

4.3 Poiseuille flow: Body forces with volume penalisation

We now test the robustness of optimal masks to forcing terms in Poiseuille flow, steady pressure driven flow between two no-slip walls. This solution tests the robustness of optimal masks to forcing terms in the problem. The penalised problem reduces to a single equation for the tangential velocity $v$ in terms of the wall-normal direction $x$,

[TABLE]

A true wall replaces damping with no-slip boundary conditions at $x=0$. In the fluid, the flow behaves parabolically, and there is an $\mathcal{O}(\varepsilon^{2})$ flow due to pressure within the solid (fig. 7). However, the same $\mathcal{O}(\varepsilon)$ displacement length occurs for conventional masks. Optimised masks eliminate the $\mathcal{O}(\varepsilon)$ displacement length, promoting the average error $\mathcal{E}_{1}$ to $\mathcal{O}(\varepsilon^{2})$, though the maximum local error $\mathcal{E}_{\infty}$ at $x=0$ remains $\mathcal{O}(\varepsilon)$ (fig. 7 ($c$)). We note that though it is possible to further reduce the fluid error outside of the mask, it is not possible to reduce the global error below $\mathcal{O}(\varepsilon^{2})$, or the local error below $\mathcal{O}(\varepsilon)$. This is because the error in the boundary region begins to dominate the error in the fluid region.

4.4 Viscous stagnation-point flow: Reynolds dependent damping

Flow about a stagnation point adds two important physical effects; flow normal and tangent to a boundary, and equation nonlinearities parameterised by the Reynolds number $\operatorname{Re}$. We can use symmetry to reduce the problem to a single third order nonlinear differential equation for the wall-normal velocity $u$ in the wall-normal coordinate $x$. Adding a finite volume penalisation boundary, we now solve for the penalised normal velocity $u$,

[TABLE]

We cannot extend our volume penalised region indefinitely, due to the linear pressure driven flow in this region, which causes the velocity at $x=0$ to increase as the true boundary is moved further away. We instead enforce no slip boundary conditions on $u$ and $u^{\prime}$ at the non-dimensional location $x=-1$. The true profile replaces damping with no-slip conditions applied at $x=0$. We solve the system numerically using Dedalus [9].

This simple system reproduces the distinction between normal and tangential velocity errors noted in section 3. The combination of no-slip boundaries and incompressibility imply zero gradient in the normal velocity at the boundary fig. 7 ($b$), ($d$). Hence the normal velocity is always $\mathcal{O}(\varepsilon^{2})$ in the boundary region, whereas the tangential velocity can be $\mathcal{O}(\varepsilon)$.

Stagnation point flow also contains nonlinearity and the Reynolds number $\operatorname{Re}$, leading to an additional damping regime. As in section 3, we see a transition from intermediate to strong damping past $\eta<\operatorname{Re}^{-1}$ in a plot of the global error $\mathcal{E}_{1}$ of $u$ as a function of $\eta$ for different $\operatorname{Re}$ in fig. 8. While unoptimised masks can achieve temporary $\mathcal{O}(\eta)$ accuracy for intermediate damping, only zero displacement length masks can achieve $\mathcal{O}(\eta)$ convergence in large damping regimes (fig. 7 ($d$)). This peculiar, if fortunate, high Reynolds number behaviour was observed experimentally by Angot et al. [3], contradicting the theoretical convergence predicted by Angot [2] ($\mathcal{O}(\varepsilon^{1/2})$) and Carbou and Fabrie [14] ($\mathcal{O}(\varepsilon)$).

5 Flow past a rotating cylinder: Moving objects in 2D

We now examine 2D incompressible flow past a rotating cylinder at moderate Reynolds number. This system includes damping, viscosity, non-linearity, curvilinear geometries, and unsteady linear and rotational acceleration of the flow and body respectively, leading to non-zero unsteady horizontal drag $F_{x}$, vertical drag $F_{y}$, and torque $T$. We simulate both a true no-slip cylinder and various volume penalised approximations using the open-source spectral PDE solver Dedalus [9] to determine the local and global error of each field as well as the physically relevant drag and torque accuracy. We compare the performance of standard unoptimised masks with optimally shifted and optimally smoothed masks, validate our prescriptions, and apply a simple Richardson extrapolation scheme to further accelerate the convergence of optimised masks to $\mathcal{O}(\varepsilon^{4})\sim\mathcal{O}(\eta^{2})$ accuracy. The configurations are tested with geometry conforming spectral elements at $\operatorname{Re}=200$, to examine the mathematical convergence in isolation of numerical error. First we examine the mathematical problem.

5.1 Mathematical formulation

We simulate the volume penalised cylinder of unit radius on an annular domain $(\theta,r)\in[0,2\pi]\times[R_{1},R_{2}]$. We do not include the center to avert the coordinate singularity. We then solve the initial value problem with zero initial conditions, and governing equations

[TABLE]

This is a slight reformulation of standard incompressible hydrodynamics in polar coordinates. We solve for the radial and polar velocity $u$ and $v$, and have introduced the variables $p$ and $q$ to write the problem in first order form (as required for Dedalus), which are related to the true pressure $P$ and vorticity $\omega_{z}$ by

[TABLE]

We also multiply out factors of $r^{-1}$ to improve the convergence of spectral discretisations of the equations. The outer boundary at $r=R_{2}$ applies uniform flow boundary conditions accelerating from rest

[TABLE]

where the time-dependent $g(t)$ utilises the compact mask function $\Gamma_{[\operatorname{erf};1]}$ to smoothly accelerate from $g(0)=0$ to $g(\Delta t)=1$ at $\Delta t=4$, with steady uniform flow boundary conditions ${g(t>\Delta t)=1}$ until the end time $t_{end}=10$. The inner coordinate singularity at $r=0$ is avoided by prescribing rotating no slip inner boundaries within the mask at $r=R_{1}=1/10$,

[TABLE]

where the cylinder angular velocity $\Omega$ is sinusoidal. The reference simulation replaces the volume penalisation with the rotating no-slip boundary conditions applied at $r=R_{1}=1$.

Force calculations — The force $F_{0}$ and torque $T_{0}$ on a no slip boundary are defined by surface integrals of the tangential and normal components of the stress tensor

[TABLE]

The penalised force $F$ and torque $T$ calculations of eq. 5 including correction terms are

[TABLE]

Azimuthal symmetry means the acceleration correction appears only for the torque.

Volume-penalty parameters— Three different types of penalty mask are chosen,

• •

$\Gamma_{0,0}$: Standard discontinuous mask $\ell=0,\delta=0$,

• •

$\Gamma_{\ell,0}$: Shifted discontinuous mask $\ell=\varepsilon,\delta=0$,

• •

$\Gamma_{0,\delta}$: Zero-shift smoothed mask $\ell=0,\delta=3.80171928\varepsilon$,

which are all defined as various limits of a smooth compactified error function

[TABLE]

We then examine these mask families for $\eta_{1}=10^{-2}$, $\eta_{2}=10^{-3}$, and $\eta_{3}=10^{-4}$.

5.2 Numerical method and parameters

The initial value problem is simulated using the Dedalus framework [9] 111Visit the Dedalus home page at http://dedalus-project.org/.. Dedalus is a general-purpose framework for solving partial differential equations using spectral methods. Dedalus allows users to construct domains from the direct product of spectral series including Fourier Series, Chebyshev polynomial series, and continuous segments of Chebyshev polynomials (a compound Chebyshev basis), and to enter systems of PDEs in plain text. The terms on the left-hand side (LHS) of the equations must be linear and are discretised into sparse and banded matrices acting on the spectral coefficients of the solution. Boundary conditions are enforced at the endpoints of the Chebyshev bases with the use of a preconditioned Chebyshev-tau method. The equations can be evolved using a range of mixed implicit-explicit timesteppers, with the LHS terms being integrated implicitly, and the right-hand side (RHS) terms integrated explicitly. The framework is written in the Python programming language but utilises compiled libraries for performance and automatically parallelises the solvers using MPI, enabling for rapid prototyping and comparisons between models as well as efficient simulations at scale.

The test and reference solutions use well-resolved cylindrical grids in Dedalus to calculate the mathematical error of the penalised problem. We use a Fourier basis of 512 modes for the azimuthal direction and Chebyshev bases in the radial direction. A single Chebyshev basis with 256 modes is used for the reference solution, over the interval $I=[1,R_{2}]$. No-slip boundary conditions are enforced at $r=1$. A compound Chebyshev basis is used for the penalised simulations, with additional segments to resolve the mask and the interior of the cylinder. We partition over the intervals $I_{1}=[R_{1},1+\ell-\delta],I_{2}=[1+\ell-\delta,1+\ell+\delta]$, and $I_{3}=[1+\ell+\delta,R_{2}]$, where $R_{1}=0.1$ and $R_{2}=10$. Each interval is discretised using 256 modes, 64 modes, and 256 modes respectively. In the singular case $\delta=0$ the middle interval is neglected.

The system is timestepped using a third order accurate backward difference scheme. Other timesteppers (including Crank-Nicolson Adams-Bashforth and Runge-Kutta schemes) performed poorly with the time varying velocity boundary conditions for $t<4$. The scheme was validated as converged by observing spectral accuracy ($\mathcal{O}(10^{-10})$ error) for the pressure equation and boundary conditions at all times during the simulation. The penalty term is timestepped explicitly (for accuracy reasons [37]) with the stability constraint suggesting our choice of timestep $5\times 10^{-5}$. This system is then timestepped to $t=10$ for all simulations. Minimal working examples of all simulations in this paper are available on GitHub at [26].

5.3 Richardson extrapolation of model error

Richardson extrapolation is a general technique to accelerate convergence of sequences with known asymptotic behaviour. For a quantity $X_{i},X_{j}$ derived from simulations with penalty parameters $\eta_{i}$ and $\eta_{j}$ we can calculate an extrapolated quantity $\overline{X}_{ij}$ as

[TABLE]

If $X_{i},X_{j}$ obey leading order behaviour proportional to $\eta$, the extrapolation $\overline{X}_{ij}$ will cancel this error. One can thus achieve $\mathcal{O}(\eta^{2})$ accuracy in practice, for the small added cost of a lower resolution simulation at larger $\eta$.

5.4 Summary of errors

Figure 9 compares the pressure $P$ and vorticity $\omega_{z}$ fields for the reference simulation and the volume-penalty simulation $(2,2)$. Snapshots are shown every two seconds until time $t=10$, and a black dot illustrates the rotation of the cylinder. The pressure color scale is calibrated to the steady state uniform flow after $t=4$, and the figures are plotted to $r=8$ to capture the vortex evolution. After $t=4$ the outer boundaries are steady and higher pressures are seen at the streamwise side of the cylinder, with low pressures around vortices and the sides of the cylinder. The pressure varies within the volume penalised cylinder, requiring the inner stress correction term in the force and torque calculations. Vortex shedding is observed immediately due to the rotation of the cylinder and the inherent instability of the von Karman vortex street.

To quantify the error fig. 10 plots the pointwise error $u-u_{0}$ in the radial velocity field at $t=10$ for the different masks ($\Gamma_{0,0},\Gamma_{\ell,0},\Gamma_{0,\delta}$) and penalty parameters $\eta=\eta_{0},\eta_{1},\eta_{2}$. The colormap range is normalised to half the maximum local error $\mathcal{E}_{\infty}$, printed below each figure. Optimal adjustments give identical spatial profiles of the error (the two right columns), with the amplitude decreasing as $\mathcal{O}(\eta)$. Standard discontinuous masks lack this property, with varying spatial pattern and reduced convergence in $\eta$. This is strong empirical evidence for the validity of the asymptotic expansions of section 3.

This motivates the extrapolated calculations in fig. 10. Extrapolations of consecutive values of $\eta$ were performed for each mask. No improvement is observed for the standard discontinuous mask, as expected from the $\mathcal{O}(\varepsilon)$ displacement length error. However extrapolation is very effective for optimised masks. Extrapolation using $\eta_{0}=10^{-2}$ and $\eta_{1}^{-3}$ performs similarly to $\eta_{2}=10^{-4}$, and extrapolation from $\eta_{1}$ and $\eta_{2}$ is almost two orders of magnitude more accurate than $\eta_{2}$ alone. This shows $\mathcal{O}(\eta^{2})$ convergence for Richardson extrapolations of optimal masks, a massive performance boost. By running a smaller second simulation in tandem with the first, it is possible to improve convergence to $\mathcal{O}(\eta^{2})\propto\mathcal{O}(\varepsilon^{4})$, implying only twice the effort to halve the error in 2D, and slightly more in 3D, vastly more efficient than the 16 and 32-fold increases for uncorrected masks. This extrapolation is performed post simulation, but could be calculated during a run. The Lyapunov time of chaotic fluid flows would require reinitializing each simulation over this time scale to prevent desynchronisation.

We then calculate the global $\mathcal{E}_{1}$ and local $\mathcal{E}_{\infty}$ error norms at the final timestep for the radial velocity $u$, the azimuthal velocity $v$, the true pressure $P$, and the true vorticity $\omega_{z}$ in fig. 11. We see the standard discontinuous mask $\Gamma_{0,0}$ (black) performs significantly worse than optimised masks, giving the predicted $\mathcal{O}(\eta^{1/2})$ convergence rate. Optimised masks instead show $\mathcal{O}(\eta)$ convergence, and extrapolated predictions achieve $\mathcal{O}(\eta^{2})$ error (based off the smallest $\eta$ used). However, for the weakest penalisation $\eta_{0}=$1\text{\times}{10}^{-2}$$, all approaches are equally accurate; In the intermediate damping regime, the dominant error is not due to the displacement length (which optimised masks lack), but time scale errors due to pressure gradients and acceleration of the penalty mask, as seen for stagnation point flow in section 4.1. Similar behaviour for the global error is observed for the azimuthal velocity $v$ and pressure $P$.

The differing behaviour of the global and local errors for each variable stems from the differentiability of the true solutions at the boundary, and the localisation of error achieved by optimal masks. The radial velocity is $C^{1}$ as it and its derivative are zero at the boundary. This is uniformly approximated by penalised solutions, giving commensurate global $\mathcal{E}_{1}$ and local $\mathcal{E}_{\infty}$ errors. However the tangential velocity $v_{0}$ and pressure $p_{0}$ are $C^{0}$ with an $\mathcal{O}(1)$ jump in the derivative at the boundary (e.g. fig. 4). This is less regular than the volume penalised solution, causing $\mathcal{O}(\varepsilon)$ disagreement in the boundary region. Calibrating the mask localises this error to the boundary, implying maximum error of $\mathcal{O}(\varepsilon)$, but global error of only $\mathcal{O}(\varepsilon^{2})$. The derivative (vorticity) is discontinuous at the boundary, causing $\mathcal{O}(1)$ disagreement in the interior of the mask near the boundary. Optimal masks cause $\mathcal{O}(1)$ disagreement in $\omega_{z}$ in an $\mathcal{O}(\varepsilon)$ boundary region and $\mathcal{O}(\varepsilon^{2})$ elsewhere, while standard masks cause $\mathcal{O}(\varepsilon)$ error everywhere.

We plot time series of the reference horizontal drag $F_{0,x}$, vertical lift $F_{0,y}$, and reference torque $T_{0}$ in fig. 12. Heightened drag is observed during the peak fluid acceleration ($t<\Delta t=4$), corresponding to the large induced pressure gradient observed at $t=2$ in fig. 12. A significant lift force is also generated during the simulation, as expected from the Magnus effect [41].

Figure 13 compares the time series of the horizontal drag error $\Delta F_{x}$, vertical drag error $\Delta F_{y}$, and torque error $\Delta T$ for each choice of mask $\Gamma_{0,0},\Gamma_{\ell,0}$ and $\Gamma_{0,\delta}$ and choice of damping $\eta=\eta_{0},\eta_{1},\eta_{2}$. All the masks improve similarly in the transition from $\eta_{1}=10^{-2}$ to $\eta_{2}=10^{-3}$, in keeping with the intermediate damping regime. However, beyond $\eta=10^{-2}$ the unshifted discontinuous mask $\Gamma_{0,0}$ transitions from $\mathcal{O}(\operatorname{Re}\varepsilon^{2})$ to $\mathcal{O}(\varepsilon)$ convergence, whereas optimised masks maintain consistent $\mathcal{O}(\eta)$ convergence. The temporal profile of the error is largely identical for the corrected masks, but inconsistent for the unshifted discontinuous mask. Applying the extrapolation procedure to the drag, lift and torque calculations, and we find similar performance improvements. No improvement occurs for the standard mask, while the adjusted masks each show $\mathcal{O}(\eta^{2})$ convergence. This is an extraordinary improvement in the accuracy of the volume-penalty method, attaining drag errors of order $10^{-5}$ using extrapolation at $\eta=10^{-4}$.

6 Conclusions

In this paper we examine and improve the volume-penalty method, a model of objects in fluids which approximates standard no-slip boundaries. A key observation is that damping induces both a penalisation time scale $\eta$, and length scale $\varepsilon$, related by the Reynolds number $\eta=\operatorname{Re}\,\varepsilon^{2}$. By using multiple-scales matched-asymptotics to second order in $\varepsilon$, and convenient signed distance boundary layer coordinates, we show that the leading order model error in the fluid is in general $\mathcal{O}(\varepsilon)$ due to the displacement length of general mask functions.

We then show this displacement error can be eliminated by correctly shifting and smoothing the mask near the boundary, promoting the accuracy to $\mathcal{O}(\eta)$. This calculation can be done once efficiently for a particular mask function using a Riccati transformation. This improvement in model error translates to an improvement in total error from first to second order numerical convergence.

The second order asymptotics also reveals two damping regimes for large $\operatorname{Re}$. Intermediate damping $1>\eta>\varepsilon>\operatorname{Re}^{-1}$ implies the ‘second order’ time scale error dominates the ‘first order’ displacement length error, while strong damping $\operatorname{Re}^{-1}>\varepsilon>\eta$ shows the displacement length error dominates and mask corrections must be used to achieve $\mathcal{O}(\eta)$ accuracy. We also derive simple force and torque calculations which are proportional to the accuracy of the flow.

We validate each of these findings in Poiseuille flow, viscous stagnation point flow, and accelerating flow past a rotating cylinder at $\operatorname{Re}=200$. We finally show that in the strong damping regime, Richardson extrapolation can also be used with optimised masks to boost model error convergence to $\mathcal{O}(\eta^{2})$.

Acknowledgments

Eric Hester acknowledges support from The University of Sydney through the Postgraduate Teaching Fellowship, and the Phillip Hofflin International Research Travel Scholarship. Geoffrey Vasil acknowledges support from the Australian Research Council, project number DE140101960. The authors would also like to acknowledge the thoughtful input from many discussions with Daniel Lecoanet, Christopher Lustri, Benjamin Favier, Laurent Duchamin, and Kai Schneider.