Volume Penalty Method

Many researchers have their preferred CFD codes for solving fluid dynamics. I happen to think that Dedalus is ideal for exploring new physics. It’s incredibly flexible, but also performant enough to run serious problems on thousands of cores. But there’s a pretty significant barrier to using spectral codes like Dedalus for simulating fluid-solid interactions, and interesting geometries more generally. Namely, spectral codes work well on tensor-product domains (cubes, spheres, sphered cubes). They don’t work well on weird geometries.

One way out is to use coordinate remappings to transform the domain into a tensor-product domain. But transforming the equations involves a lot of differential geometry, and poses all sorts of new stability problems numerically. Fortunately, there is an easier way.

Turn boundary conditions into equation source terms!

This simple idea underlies many numerical methods, like phase-field models, diffuse-interface methods, diffuse-domain methods, fictitious-domain methods, immersed-boundary methods, and of course, the volume-penalty method. In the case of the volume-penalty method, we are approximating no-slip solid boundaries in fluid dynamics with rapid linear damping in the “fictitious” solid interior. Mathematically, this looks like dropping the annoying solid boundary conditions, and replacing them with simple source terms in the Navier-Stokes equations: $$ \begin{align} \partial_t u + u\cdot \nabla u + \nabla p -\nu \nabla^2 u &= 0, \qquad u = 0 \text{ on } \partial \Omega_s,\\ \text{turns in to}&\\ \partial_t u + u\cdot \nabla u + \nabla p -\nu \nabla^2 u &= -\frac{\Gamma}{\eta} u \end{align} $$ where $\Gamma$ is approximately an indicator function over the solid $\Omega_s$, and $\eta\ll 1$ is a small damping timescale.

This simple method is how I simulated the flowing dye around an ellipse in the header video. The question is, how accurate is this method, and can it be improved? That’s what my first paper is all about.