# A volume penalization method for incompressible flows and scalar advection-diffusion with moving obstacles

Benjamin Kadoch, Dmitry Kolomenskiy, Philippe Angot, Kai Schneider

► To cite this version:

Benjamin Kadoch, Dmitry Kolomenskiy, Philippe Angot, Kai Schneider. A volume penalization method for incompressible flows and scalar advection-diffusion with moving obstacles. Journal of Computational Physics, 2012, 231 (12), pp.4365-4383. 10.1016/j.jcp.2012.01.036. hal-01032208

HAL Id: hal-01032208

https://hal.science/hal-01032208

Submitted on 15 Mar 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# A volume penalization method for incompressible flows and scalar advection–diffusion with moving obstacles

Benjamin Kadoch<sup>a,b,*</sup>, Dmitry Kolomenskiy<sup>a,c</sup>, Philippe Angot<sup>e,d</sup>, Kai Schneider<sup>a,d</sup>

<sup>a</sup> M2P2–UMR 6181 CNRS, Aix-Marseille Université, Ecole Centrale Marseille, Marseille, France

<sup>b</sup> IUSTI–UMR 6595 CNRS, Aix-Marseille Université, Marseille, France

<sup>c</sup> CERFACS, Toulouse, France

<sup>d</sup> Centre de Mathématiques et d'Informatique, Aix-Marseille Université, Marseille, France

<sup>e</sup> LATP–UMR 6632 CNRS, Aix-Marseille Université, Marseille, France

A volume penalization method for imposing homogeneous Neumann boundary conditions in advection–diffusion equations is presented. Thus complex geometries which even may vary in time can be treated efficiently using discretizations on a Cartesian grid. A mathematical analysis of the method is conducted first for the one-dimensional heat equation which yields estimates of the penalization error. The results are then confirmed numeri-cally in one and two space dimensions. Simulations of two-dimensional incompressible flows with passive scalars using a classical Fourier pseudo-spectral method validate the approach for moving obstacles. The potential of the method for real world applications is illustrated by simulating a simplified dynamical mixer where for the fluid flow and the sca-lar transport no-slip and no-flux boundary conditions are imposed, respectively.

## 1. Introduction

Computational fluid dynamics (CFD) in complex geometries which may vary in time, a problem typically encountered in fluid-structure interaction problems, is still a major challenge and requires advanced numerical techniques. Different approaches have been developed so far and an exhaustive overview is beyond the scope of the paper. Body fitted coordinate systems, see, e.g., [1], yield a well adapted discretization for a given geometry. For time varying geometries the grid generation becomes however even more complex and for instance elliptic grid generation techniques [2,3] are necessary which further increase the computational effort.

During the last decades immersed boundary methods gained ground and became an attractive alternative. The idea can be traced back to Courant (1943) [4] in the context of constrained optimization to obtain problems free of constraints. There exists a large variety of immersed boundary methods, for example Lagrangian multipliers [5], level-set methods [6], fictitious domain approaches and surface [7] and volume penalization approaches [8]. For reviews we refer the reader to [9,10].

A common feature of immersed boundary methods is that the Navier–Stokes equations can be discretized on a classical grid, e.g., a Cartesian one in a simple domain, for instance, a square, for which efficient solvers are available. The complex geometry is then immersed into the larger computational domain and the boundary conditions are appended by adding terms to the equations. The imposed boundary conditions are classically of Dirichlet type corresponding to no-slip walls for fluid flow. For Neumann boundary conditions, which correspond to no-flux conditions only some works can be found

* Corresponding author at: IUSTI–UMR 6595 CNRS, Aix-Marseille Université, Marseille, France.

E-mail addresses: benjamin.kadoch@polytech.univ-mrs.fr (B. Kadoch), dkolom@L3m.univ-mrs.fr (D. Kolomenskiy), angot@cmi.univ-mrs.fr (P. Angot), kschneid@cmi.univ-mrs.fr (K. Schneider).
in the literature, for example [11,12]. A drawback of most immersed boundary methods is their limitation to low order approximation, i.e., first or second order, and their computational stiffness leading to ill-conditioned linear systems or prohibitive small time steps.

In the current work we focus on the volume penalization approach which is physically motivated by modeling solid walls as porous media whose permeability tends to zero [13] and mathematically justified in [14,15]. Therewith no-slip boundary conditions have been applied in many simulations [8,16,17], including benchmarking of dipole-wall interactions [18,19] and turbulent flows in confined geometries [20,21]. The volume penalization has been generalized to impose in addition to Dirichlet [8] also Neumann or Robin boundary conditions [22]. However, this generalization has been applied in the finite element or finite volume context but is not directly applicable to pseudo-spectral methods.

The aim of the present paper is to propose a method well adapted to study transport and mixing in complex geometries. The novelty is to include into a classical Fourier pseudo-spectral code, homogeneous Neumann boundary conditions for advection-diffusion of passive scalars combined with the Navier-Stokes equations describing the fluid flow. Another challenge is to take into account time-dependent geometries with homogeneous Neumann boundary conditions using the same approach as in [23] which generalizes the volume penalization method for moving obstacles with Dirichlet boundary conditions. Since spectral methods are available for direct numerical simulation of turbulence in periodic domains, the method proposed here could be implemented easily into existing codes. A mathematical analysis will be conducted in order to validate the penalization approach for Neumann boundary conditions which is confirmed by numerical simulations. An example for more complex flows in a dynamical mixer illustrates the potential of the method.

Studying scalar mixing is motivated by many applications, such as chemical reactors, ventilation systems or atmospheric pollutant mixing. Numerical simulation is a powerful tool to characterize these industrial devices such as mixers, see, e.g., [24]. Numerous studies have been carried out also in the context of homogeneous isotropic turbulence, see, e.g., [25,26] or in presence of walls [27,28]. However, almost all the flows in nature and in industrial devices are wall bounded, so we need to include solid boundaries to be able to simulate these configurations. In recent work two of the authors studied the influence of wall-bounded domains on the Lagrangian statistics of particles [21] and they showed that in particular high order statistics differ significantly from those in periodic domains. This motivates the study of passive scalars in the presence of walls.

The manuscript is decomposed as follows. In Section 2 the mathematical method will be described presenting the governing equations. A mathematical analysis of the penalization technique is then given for simplified examples. Section 3 shows numerical simulations in one and two space dimensions which are in agreement and thus validate the theoretical results. Then in Section 4 the penalization method will be extended to moving obstacles and illustrations for different two-dimensional flows with passive scalars will be presented. Finally, some conclusions and perspectives will be exposed.

## 2. Mathematical method

The computational domain $\Omega$ contains the fluid domain $\Omega_f$ and the solid domain $\Omega_s$, as illustrated in Fig. 1. The volume penalization method consists in considering the whole domain $\Omega$ and to modify the governing equations in order to insure the desired boundary conditions at the interface $\Sigma$ between the fluid and solid domains. In the case of spectral methods, the computational domain $\Omega$ is considered to be periodic.

The fluid flow is governed by the incompressible Navier-Stokes equations, which we recast here in the form of the penalized momentum equation:

$$ \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} - \nu \nabla^2 \mathbf{u} + \frac{1}{\rho} \nabla p = - \frac{1}{\eta_u} (\chi \mathbf{u}) \quad \text{for } x \in \Omega, t > 0 \quad (1) $$

<figure><img src="image_2.png" /><figcaption>Fig. 1. The computational domain $\Omega$ is the union of the fluid domain $\Omega_f$ and the solid domain $\Omega_s$. The interface between these two domains is $\Sigma = \partial\Omega_f \cap \partial\Omega_s$ and $\Gamma = \partial\Omega$ is the boundary of $\Omega$.</figcaption></figure>
where the velocity **u** is divergence-free, $\nabla \cdot \mathbf{u} = 0$, $p$ is the pressure, $\rho$ is the fluid density assumed to be constant with $\rho = 1$ and $v$ the kinematic viscosity. The right-hand side of Eq. (1) corresponds to the penalization term, which is used to impose the no-slip boundary condition, $\mathbf{u}|_{\Sigma} = 0$. The mask function $\chi$ equals 1 inside the solid $\Omega_s$ and 0 inside the fluid domain $\Omega_f$, where the original (non penalized) momentum equation is recovered. The permeability $\eta_u$ is required to be sufficiently small for a given $v$ in order to minimize the modeling error [8,17], since the physical idea behind the volume penalization method is to consider solids as porous media with vanishing permeability, such that the velocity of the surrounding fluid vanishes at the interface. Angot et al. [8] proved rigorously the convergence of the velocity of the penalized equation in the limit of vanishing $\eta_u$ to the solution of the velocity of the Navier-Stokes with no-slip boundary conditions, with a global convergence rate as $O(\eta_u^{1/4})$. Carbou and Fabrie [15] improved this estimate by using BKW analysis of the boundary layer for smooth solutions. In the $H^1$-norm, the convergence rate varies as $O(\eta_u^{1/4})$ inside the solid and as $O(\eta_u^{1/2})$ in the fluid domain, whereas the error $L^2$-norm varies like $O(\eta_u^{3/4})$ inside the solid and $O(\eta_u^{1/2})$ in the fluid domain.

To study mixing, we extend the penalization technique to the advection-diffusion equation of a passive scalar:

$$ \frac{\partial \theta}{\partial t} + ((1-\chi)\mathbf{u}) \cdot \nabla \theta = \nabla \cdot ([\kappa(1-\chi) + \eta_\theta \chi] \nabla \theta), \quad (2) $$

where $\theta$ is the passive scalar, which can represent for example the concentration or the temperature, $\kappa$ is the diffusivity, and $\eta_\theta$ is the penalization parameter. The right-hand side of Eq. (2) corresponds to the penalized diffusion term [22], which takes into account the no-flux boundary condition at the wall $\nabla\theta \cdot \mathbf{n}|_{\partial\Omega} = 0$ where $\mathbf{n}$ is the normal vector to the wall. Zero flux through the fluid-solid interface can be achieved by imposing vanishing diffusivity inside the solid domain. In the periodic case without boundaries where the penalization term is absent this term simply equals $\kappa \nabla^2\theta$.

In the following analysis we neglect the advective term in Eq. (2), that is we assume that $\mathbf{u} \equiv 0$.

## 2.1. One-dimensional diffusion equation

In this part we study the simple test case, in which the penalized problem has an analytical solution. We consider a one-dimensional diffusion equation on the positive real axis with the Neumann boundary condition imposed at the origin,

$$ \frac{\partial \theta}{\partial t} - \frac{\partial^2 \theta}{\partial x^2} = 0, \quad x > 0, \quad t > 0 \quad (3) $$

completed with the initial and boundary conditions

$$ \begin{cases} \theta|_{t=0} = \cos x, \\ \frac{\partial \theta}{\partial x}|_{x=0} = 0, \end{cases} \quad (4) $$

and where $\theta$ is supposed to be a bounded function. The exact solution of this problem is

$$ \theta(x, t) = e^{-t} \cos x. \quad (5) $$

Now we model the Neumann boundary condition in Eq. (3) using the volume penalization method. For convenience of notations, in this part we drop the subscript of $\eta_\theta$ and simply denote it by $\eta$. The spatial domain is extended to the whole real axis, and for negative $x$ the solution is governed by the diffusion equation with diffusivity $\eta$, i.e., the penalization parameter, which is small. We shall call $x > 0$ the 'fluid domain', and $x < 0$ the 'solid domain'. Continuity of the solution and of the flux is imposed at the origin $x = 0$ for $t > 0$. The initial condition in the solid domain is zero. Hence the penalized problem reads,

$$ \left\{ \begin{aligned} \frac{\partial \theta_{\eta}}{\partial t} - \frac{\partial^2 \theta_{\eta}}{\partial x^2} &= 0, & x > 0, & t > 0, \\ \frac{\partial \theta_{\eta}}{\partial x} - \eta \frac{\partial^2 \theta_{\eta}}{\partial x^2} &= 0, & x < 0, & t > 0, \\ \theta_{\eta}|_{t=0} &= \cos x, & x > 0, \\ \theta_{\eta}|_{t=0} &= 0, & x \le 0, \\ \lim_{x \to 0^{-}} \theta_{\eta} &= \lim_{x \to 0^{+}} \theta_{\eta}, \\ \lim_{x \to 0^{-}} \eta \frac{\partial \theta_{\eta}}{\partial x} &= \lim_{x \to 0^{+}} \frac{\partial \theta_{\eta}}{\partial x}. \end{aligned} \right. \quad (6) $$

This problem can be solved using the same approach as previously applied for Dirichlet boundary conditions [16,29]. Performing the Laplace transform in time, $v_\eta(\alpha, x) = \int_0^\infty e^{-\alpha x} \theta_\eta(x, t) dt$, $\alpha \in \mathbb{C}$, we obtain the following set of ordinary differential equations:

$$ \left\{ \begin{aligned} \alpha v_{\eta} - \cos x - v_{\eta}'' &= 0, & x > 0, \\ \alpha v_{\eta} - \eta v_{\eta}'' &= 0, & x < 0, \end{aligned} \right. \quad (7) $$
where primes stand for derivatives with respect to x. These equations can readily be solved, and integration constants are determined by requiring $v_η$ to be bounded at infinity and using the continuity conditions. Finally, after performing the inverse Laplace transform, we obtain

$$ \left\{ \begin{aligned} \theta_{\eta}(x, t) &= e^{-t} \cos x - \frac{1}{2\sqrt{\pi}} \frac{\sqrt{\eta}}{1 + \sqrt{\eta}} \int_0^t \frac{e^{\tau-t} x e^{-x^2/4\tau}}{\tau^{3/2}} d\tau, & x > 0, \\ \theta_{\eta}(x, t) &= -\frac{1}{2\sqrt{\pi}} \frac{1}{1 + \sqrt{\eta}} \int_0^t \frac{e^{\tau-t} x e^{-x^2/4\eta\tau}}{\sqrt{\eta}\tau^{3/2}} d\tau, & x < 0. \end{aligned} \right. \qquad (8) $$

The cosine in the right-hand side of the first of these equations is recognizable as Eq. (5), the exact solution of Eq. (3). The remaining term gives the penalization error in the fluid domain $x > 0$, $e_f = \theta_η - \theta$. This error is shown in Fig. 2 (left) as a function of x at three different time instants and for two different values of $\eta$. The plot on the right shows its derivative with respect to x. Both are decaying with $x \to \infty$ and with $t \to \infty$. Smaller values of $\eta$ result in smaller errors. The error in the solid domain, $x < 0$, is identical to the solution of the penalized problem, $e_s = \theta_η$. It is shown in the left panel of Fig. 3, while its derivative is depicted in the right panel. This figure reveals a thin boundary layer inside the solid, similar to that observed in [16] for the Dirichlet boundary condition. While $e_s$ is large (compared to $\max|e_f|$) and $\partial e_s / \partial x$ is even larger, the thickness of this layer diminishes with decreasing $\eta$.

The spatial $L^2$-norm of the error in the fluid domain and the norm of the derivative equal, respectively,

$$ \|e_f\|_{L^2(\Omega_f)} = \frac{\sqrt{\eta}}{1 + \sqrt{\eta}} F_1(t), \quad \|\partial e_f / \partial x\|_{L^2(\Omega_f)} = \frac{\sqrt{\eta}}{1 + \sqrt{\eta}} F_2(t). \qquad (9) $$

Correspondingly, in the solid domain we have,

$$ \|e_s\|_{L^2(\Omega_s)} = \frac{\eta^{1/4}}{1 + \sqrt{\eta}} F_1(t), \quad \|\partial e_s / \partial x\|_{L^2(\Omega_s)} = \frac{\eta^{-1/4}}{1 + \sqrt{\eta}} F_2(t). \qquad (10) $$

with the time dependent functions,

$$ F_1(t) = \left\{ \int_0^\infty \left( \frac{\partial f(y, t)}{\partial y} \right)^2 dy \right\}^{1/2}, \quad F_2(t) = \left\{ \int_0^\infty \left( \frac{\partial^2 f(y, t)}{\partial y^2} \right)^2 dy \right\}^{1/2}, \qquad (11) $$

and where

$$ f(y, t) = \int_{0}^{t} \frac{e^{\tau-t} e^{-y^2/4\tau}}{\sqrt{\pi \tau}} d\tau. \qquad (12) $$

The graphs of the time dependent functions $F_1$ and $F_2$ are depicted in Fig. 4. Thus, in the fluid domain, $\|e_f\|_{H^1}^2 = \|e_f\|_{L^2}^2 + \|\partial e_f / \partial x\|_{L^2}^2$ converges to zero with $\eta \to 0$, and it also decays when time becomes large. However, when $t \to 0$, the error in the derivative is diverging. This is a consequence of the discontinuity in the initial condition. Care is needed, however, when choosing a continuous initial condition, since a non-zero $\theta_η|_{t=0}$ in the solid domain is acting as a source of error for the fluid domain. In the solid domain, the norm $\|e_s\|_{H^1}^2 = \|e_s\|_{L^2}^2 + \|\partial e_s / \partial x\|_{L^2}^2$ diverges with $\eta \to 0$, whereas $\|e_s\|_{L^2} = O(\eta^{1/4})$.

The cosine initial condition in Eq. (6) verifies the Neumann condition at the boundary $x = 0$. It is also important to consider a situation when the initial condition fails to satisfy the boundary condition, as this might often be the case in practical numerical simulations. Starting from the initial condition

<figure><img src="image_3.png" /><figcaption>Fig. 2. Penalization error $e_f = \theta_η - \theta$ in the fluid domain (left); derivative of the penalization error (right).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Penalization error $e_s = \theta_\eta$ in the solid domain (left); derivative of the penalization error (right).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Time evolution of the norm of the functions $f$, and its first and second derivatives, $\partial f/\partial y$, $\partial^2 f/\partial y^2$, respectively, defined in Eqs. (17) and (11).</figcaption></figure>

$$ \begin{cases} \theta_{\eta}|_{t=0} = \sin x, & x > 0, \\ \theta_{\eta}|_{t=0} = 0, & x \le 0, \end{cases} \qquad (13) $$

we obtain the following exact solution of the penalized equation at $t > 0$:

$$ \left\{ \begin{aligned} \theta_{\eta}(x, t) &= e^{-t} \sin x + \frac{1}{1 + \sqrt{\eta}} \int_{0}^{t} \frac{e^{\tau-t} e^{-x^2/4\tau}}{\sqrt{\pi \tau}} d\tau, & x > 0, \\ \theta_{\eta}(x, t) &= \frac{1}{1 + \sqrt{\eta}} \int_{0}^{t} \frac{e^{\tau-t} e^{-x^2/4\tau}}{\sqrt{\pi \tau}} d\tau, & x < 0. \end{aligned} \right. \qquad (14) $$

The norm of the error in the fluid domain is

$$ \|e_f\|_{L^2(\Omega_f)} = \frac{\sqrt{\eta}}{1 + \sqrt{\eta}} F_0(t), \quad \|\partial e_f / \partial x\|_{L^2(\Omega_f)} = \frac{\sqrt{\eta}}{1 + \sqrt{\eta}} F_1(t), \qquad (15) $$
and in the solid domain it equals

$$ \|e_s\|_{L^2(\Omega_s)} = \frac{\eta^{1/4}}{1 + \sqrt{\eta}} F_0(t), \quad \|\partial e_s / \partial x\|_{L^2(\Omega_s)} = \frac{\eta^{-1/4}}{1 + \sqrt{\eta}} F_1(t), \qquad (16) $$

where

$$ F_0(t) = \left\{ \int_0^\infty (f(y, t))^2 dy \right\}^{1/2}. \qquad (17) $$

The time evolution of $F_0$ is shown in Fig. 4. Note that the $H^1$-norm now remains finite for $t \to 0$, since the initial condition is continuous at $x = 0$. The rate of convergence with respect to $\eta$ is the same as in the previous case and is therefore independent of the initial condition.

## 2.2. Theorem on convergence and error estimates

The above calculations of the exact penalization error for a simplified one-dimensional diffusion equation can be generalized. In the case of fixed rigid obstacles, a short sketch of this rather technical proof is given in the appendix. For details we refer to [30].

**Theorem 2.1 (Convergence and error estimates).** With usual regularity assumptions for the data, the sequence $(0_\eta)_f = (1 - \chi)0_\eta$, solution of the penalized unsteady advection-diffusion problem, strongly converges in $L^\infty(0,T; L^2(\Omega_f))$ when $\eta \to 0$ to the solution $0_f$ of the original problem in the fluid domain $\Omega_f$ with $0_s = 0$ in the solid domain $\Omega_s$.

Moreover, the error $(0_\eta)_f - 0$ in the fluid domain for the $L^\infty(0,T; L^2(\Omega_f))$ and $L^2(0,T; H_0^1(\Omega_f))$ norms and the error on the Neumann boundary condition on $\Sigma$ for a suitable norm both converge like $O(\sqrt{\eta})$.

This theorem justifies thus the use of the volume penalization method to impose Neumann boundary conditions for advection-diffusion problems for any fixed geometry. In the next sections, the extension of this method to higher dimensions is described. Numerical validations for fixed and moving obstacles in two space dimensions are also presented.

# 3. Numerical validation

## 3.1. Numerical discretization

The fluid domain is immersed in a periodic domain of size $2\pi$. A classical Fourier pseudo-spectral discretization on a Cartesian grid is used in space. The time discretization is semi-implicit and of second order using Adams-Bashforth for the non-linear term and the penalization term, and Euler backwards for the linear diffusion term [17]. The stability limit of the time discretized Eq. (1) is imposed by the Courant-Friedrichs-Lewy (CFL) condition $\Delta t \le C\Delta x/U_{max}$, where $C < 1$ is a constant and $U_{max}$ the maximum velocity on the grid. As the penalization term is treated explicitly, the time step must also be smaller than the penalization parameter $\Delta t < \eta_u$.

Applying the same time scheme to Eq. (2) with the diffusive term being treated explicitly, the stability limit is given by $\Delta t \le C_2 \min((\Delta x)^2/\kappa, (\Delta x)^2/|\eta_0|)$ where $C_2 < 1$ is a constant. As $\kappa \gg \eta_0$, there is no stability limit imposed by $\eta_0$ which can thus be chosen as small as desired. Therefore, the choice of $\Delta t$ is only limited by the penalization parameter of the velocity $\eta_u$ and the above CFL condition.

## 3.2. One-dimensional diffusion equation

To get further insight into the numerical convergence properties of the volume penalization method for imposing the Neumann boundary condition, we first consider the one-dimensional diffusion equation. In this part, for convenience we denote again $\eta_0$ simply by $\eta$. The equation under consideration is

$$ \frac{\partial \theta}{\partial t} - \kappa \frac{\partial^2 \theta}{\partial x^2} = 0. \qquad (18) $$

The corresponding penalized equation is

$$ \frac{\partial \theta_{\eta}}{\partial t} - \frac{\partial}{\partial x} \left[ (\eta \chi + \kappa(1-\chi)) \frac{\partial \theta_{\eta}}{\partial x} \right] = 0, \qquad (19) $$

where $\chi$ is the mask function:

$$ \chi = \begin{cases} 1, & x \in [-2, -1[ \cup [1, 2] \\ 0, & x \in [-1, 1] \end{cases} \qquad (20) $$

The solid and fluid domains correspond thus to the domain where $\chi = 1$ and $\chi = 0$, respectively.

The computational periodic domain is $\Omega = \Omega_s \cup \Omega_f = [-2,2]$ where $\Omega_f = [-1,1]$ is the fluid domain and $\Omega_s$ is the solid domain. Homogeneous Neumann boundary conditions are imposed at $x = -1$ and $x = 1$, i.e, $\partial\theta_\eta/\partial x = 0$, for $x = \pm 1$. The initial condition, which is continuous and does not vanish inside the solid domain, is
$$ \theta(x, t=0) = \cos 4\pi x + \cos \pi x, \quad \text{for } x \in \Omega \qquad (21) $$

and it is advanced in time up to $t_{max} = 0.5$. This is done using exponential propagation (exact integration of Eq. (18) with respect to time):

$$ \theta_{\eta}(x, t_{\max}) = e^{t_{\max} \frac{\partial}{\partial x} [\eta \chi + \kappa(1-\chi)]} \theta(x, 0), \qquad (22) $$

Therefore there is no time discretization error in the numerical solution. However there is a spatial discretization error, because the spatial derivative $\frac{\partial}{\partial x}$ is approximated with the collocation differentiation matrix. We note that the exact time integration corresponds to computing the exponential of a matrix. There is also a penalization error, since $\eta$ is finite.

A series of numerical simulations is performed for $N = 2^3...2^{13}$ grid points and $\eta = 10^{-1}...10^{-8}$. The results are then compared to the exact solution of Eq. (18) at the corresponding grid points

$$ \theta(x, t) = e^{-16\pi^2 kt} \cos 4\pi x + e^{-\pi^2 kt} \cos \pi x, \qquad (23) $$

which lie inside the fluid domain $\Omega_{int} = \{x:|x|<1\}$. Subsequently, the $L^\infty$ and $L^2$-norms of the error between $\theta_{\eta}$ and $\theta$ are computed. Fig. 5 shows that the penalization error is decaying as $\eta^{1/2}$. The spatial discretization error is decaying as $N^{-1/2}$ (see Fig. 6). This is due to the low regularity of the exact solution of the penalized equation near the wall. There the solution is continuous, but the first derivative has a jump. Hence the convergence rate of the spectral method is limited. As an example, Fig. 7 shows a zoom of the exact solution Eq. (8) near the boundary at two time instants and for two values of $\eta_0$.

Fig. 8 shows the numerical solution calculated with $N = 2^{12}$ grid points for two values of $\eta$. For $\eta = 10^{-3}$, the penalization boundary layer is thick, so that there are enough grid points inside to ensure a good approximation of the penalized solution. However, the penalization error is also large. For $\eta = 10^{-6}$, there are only three grid points in the penalization boundary layer, and Gibbs oscillations are present, similar to what is observed for a discontinuous function. The $L^\infty$-error is therefore large (just slightly smaller than for $\eta = 10^{-3}$). But on average, this oscillating solution matches the exact one (continuous line), therefore the $L^2$-norm is small (much smaller than for $\eta = 10^{-3}$).

To quantify the influence of the initial condition, we performed a series of new computations using an initial condition which vanishes inside the solid domain $\Omega_s$ and inside the fluid domain, we take $\theta(x, t=0) = \cos 4\pi x + \cos \pi x$ for $x \in \Omega_f$. The results are shown in Figs. 5 and 6. The blue lines correspond to the modified initial condition (zero inside the solid) and the black lines correspond to the original initial condition (non-zero inside the solid). Actually, the discrepancy is now larger for the largest value of the penalization parameter $\eta$, because the new initial condition has a discontinuity that activates the diffusion. However, the discrepancy becomes invisible when $\eta$ is small.

An additional series of computations has been performed by changing the position of the collocation points. The latter have been shifted by $\Delta x/2$, where $\Delta x$ is the grid step size, so that the boundary lies exactly between two collocation points. The errors between the new numerical solution and the exact solution have been calculated. They are shown in Figs. 5 and 6

<figure><img src="image_6.png" /><figcaption>Fig. 5. (Color online) Decay of the approximation errors in the numerical solution of the penalized one-dimensional diffusion equation at $t = 0.5$ as a function of $\eta_0$ for the original case (black), for the collocation points shifted by a half step size (red) and for the initial condition which vanishes inside the solid domain (light blue). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Decay of the approximation error in the numerical solution of the penalized one-dimensional diffusion equation at t = 0.5 with respect to the number of Fourier modes for different but fixed η for the original case (black), for the collocation points shifted by a half step size (red) and for the initial condition which vanishes inside the solid domain (light blue). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Exact solution of the penalized one-dimensional diffusion equation Eq. (8), near the boundary x = 0.</figcaption></figure>

(red lines). The differences with respect to the original results are minor. Hence we can deduce that the relative position of the collocation points has no significant influence on the numerical precision.

We can conclude that the volume penalization method applied to the viscous term with the Neumann boundary condition thus converges and yields for sufficiently small choices of the penalization parameter accurate results.

### 3.3. Two-dimensional diffusion equation

In this part, we test the volume penalization in a two-dimensional configuration. For that, we construct a function which satisfies the two-dimensional diffusion equation with Neumann boundary conditions, and which can be used to validate the penalization method. The following problem is considered:

$$ \frac{\partial \theta}{\partial t} - \kappa \left( \frac{\partial^2 \theta}{\partial x_1^2} + \frac{\partial^2 \theta}{\partial x_2^2} \right) = 0, \quad (x_1, x_2) \in \Omega_f, \quad t \in \mathbb{R}^+, \quad \Omega_f = [-1, 1]^2, \qquad (24) $$
<figure><img src="image_9.png" /><figcaption>Fig. 8. A zoom on the solution of the one-dimensional diffusion equation near the boundary x = 1. The dashed and dashed-dotted lines show the numerical solution of the penalized equation with two different values of the penalization parameter η, as described in Section 3.2. The black line shows the exact solution of the equation with the Neumann boundary condition.</figcaption></figure>

with boundary conditions

$$ \left. \frac{\partial \theta}{\partial x_1} \right|_{x_1 = \pm 1} = 0, \quad \left. \frac{\partial \theta}{\partial x_2} \right|_{x_2 = \pm 1} = 0, \qquad (25) $$

and initial condition

$$ \theta_0(x_1, x_2) = (\cos 4\pi x_1 + \cos \pi x_1) \cos \pi x_2, \quad \text{for } (x_1, x_2) \in \Omega_f. \qquad (26) $$

The initial condition $\theta_0$ can be extended periodically, and, due to symmetry, the corresponding solution $\theta$ still satisfies the Neumann boundary condition. The exact solution is calculated as a convolution of the heat kernel with the initial condition. The Fourier transform of Eq. (24) yields

$$ \hat{\theta}(k_1, k_2, t) = e^{-vt(k_1^2+k_2^2)} \hat{\theta}_0(k_1, k_2). \qquad (27) $$

For $\theta_0$ in Eq. (26) which has been extended periodically, we thus get

$$ \begin{aligned} \hat{\theta}_0(k_1, k_2) &= \frac{1}{4\pi^2} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \theta_0(x_1, x_2) e^{-i(k_1x_1+k_2x_2)} dx_1 dx_2 \\ &= \left( \frac{\delta(k_1 - 4\pi) + \delta(k_1 + 4\pi)}{2} + \frac{\delta(k_1 - \pi) + \delta(k_1 + \pi)}{2} \right) \frac{\delta(k_2 - \pi) + \delta(k_2 + \pi)}{2}. \end{aligned} $$

Upon performing the inverse Fourier transform of Eq. (27),

$$ \theta(x_1, x_2, t) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \hat{\theta}(k_1, k_2, t) e^{i(k_1x_1+k_2x_2)} dk_1 dk_2, $$

we obtain the solution of Eqs. (24)–(26) in physical space,

$$ \theta(x_1, x_2, t) = e^{-2\pi^2 \kappa t} \cos \pi x_2 (\cos \pi x_1 + e^{-15\pi^2 \kappa t} \cos 4\pi x_1). \qquad (28) $$

In order to solve Eqs. (24)–(26) numerically using the volume penalization method, we consider an extended spatial domain $\Omega = [-1, 1]^2$, in which $\Omega_f$ is included. The mask function which defines the desired position of the boundary is

$$ \chi(x_1, x_2) = \begin{cases} 0 & \text{for } (x_1, x_2) \in [-1, 1]^2 \\ 1 & \text{elsewhere.} \end{cases} \qquad (29) $$

In the following, the solid and fluid domains correspond thus to the domains where $\chi = 1$ and $\chi = 0$, respectively.

As initial condition we use Eq. (26) evaluated in the larger domain $\Omega$. The diffusivity equals $\kappa = 10^{-2}$ and the penalization parameter is $\eta = 10^{-8}$. The discretization parameters are $N_1 \times N_2 = 512^2$ and $\Delta t = 10^{-3}$.

Fig. 9 displays the numerical solution at $t = 1$, whereas Fig. 10 shows a comparison of its one-dimensional cuts along the x- and y-axes together with the exact solution given by Eq. (28).

Moreover, the decay of the error between the exact solution and the numerical solution of the penalized equation as a function of $\eta$ is shown in Fig. 11 (black curves) for two different resolutions, $N_1 \times N_2 = 512^2$ and $1024^2$. One observes that the error decays with decreasing $\eta$. However, for sufficiently small $\eta$ the error saturates which is reflected in the plateau.
<figure><img src="image_10.png" /><figcaption>Fig. 9. (Color online) Numerical solution of the penalized two-dimensional diffusion equation at $t = 1$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Two-dimensional diffusion equation: comparison between the numerical (full line) and exact (dashed line) solutions at $t = 1$. Shown are one-dimensional cuts, one along the x-axis at $y = -0.0021$ (left) and one along the y-axis at $x = -0.0021$ (right).</figcaption></figure>

Increasing the resolution leads to smaller errors, and the saturation occurs for smaller values of $\eta$. The values of the plateau correspond to the spatial discretization error which dominates the penalization error for sufficiently small $\eta$.

The mask function is also modified by translating the discontinuity on the left by one space step to the left and the one on the right by one space step to the right, in order to quantify the influence of the relative position of collocation points. Fig. 11 (red curves) shows the computations with the new shifted mask function. The decay of the $L^2$ and $L^\infty$-errors for decreasing $\eta$ with the new shifted mask function exhibits a similar behavior as the one for the original computation and a vertical shift of the errors for the new mask functions is observed. This can be explained by the increase of the distance between the position of the collocation point defined by the position of the mask function and the real position of the interface between the fluid and solid domains $|x| = 1$. In the following, only the original mask function will be used.

Finally, it can be observed that the error decays as $N^{-3/4}$ in both the $L^2$ and the $L^\infty$-norm, which is shown in Fig. 12 for $\eta = 10^{-8}$.

Additionally we performed computations where the initial condition is changed by imposing zero inside the solid domain $\theta(x_1, x_2, t = 0) = (\cos 4\pi x_1 + \cos\pi x_2)\cos\pi x_2$, for $(x_1, x_2) \in \Omega_f$ and $0$ for $(x_1, x_2) \in \Omega_s$ (results not presented). It can be shown that if we compare with the computations done using the initial condition in Eq. (26) (black curves in Fig. 11), we observe qualitatively the same behavior as in the one-dimensional case: the errors ($L^2$ and $L^\infty$) decay with decreasing $\eta$, and then saturate for the same sufficiently small value of the penalization parameter.
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Decay of the approximation error in the numerical solution of the penalized two-dimensional diffusion equation at $t = 1$ for the original case (black) and for the new shifted mask function (red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Errors in $L^2$ and $L^\infty$-norm versus the number of Fourier modes for fixed $\eta$ at $t = 1$.</figcaption></figure>

The above computations yield thus a validation of the volume penalization method for the two-dimensional diffusion equation with Neumann boundary conditions. The precision depends on the penalization parameter which has to be chosen sufficiently small. We can note that the behavior of the solution in the solid domain is unpredictable and this is due to numerical reasons. Indeed, the solution adapts itself in the solid domain in order to guarantee the Neumann boundary condition at the boundary between the solid and fluid domain.

# 4. Moving obstacles

In this section, we extend the volume penalization method in order to compute flows containing moving obstacles. In [23], this was done for Dirichlet boundary conditions. The aim of this part is to adapt the advection-diffusion equation to
moving solids with zero flux boundary conditions at the walls. First, we validate the method using a simple test case and then we present some results for a flow in a dynamical mixer.

## 4.1. Validation

### 4.1.1. Numerical method

A two-dimensional incompressible turbulent flow with unit density is considered. The Navier-Stokes equations written in vorticity-velocity formulation are modified by a volume penalization method which takes into account the no-slip boundary conditions and the moving solid:

$$ \frac{\partial \omega}{\partial t} + (\mathbf{u} + \mathbf{U}_0) \cdot \nabla \omega - \nu \nabla^2 \omega = - \frac{1}{\eta_u} \nabla \times [\chi(\mathbf{u} + \mathbf{U}_0 - \mathbf{U}_s)], \quad (30) $$

where $\mathbf{u} = (u_1, u_2)$ is the velocity, $\omega = \nabla \times \mathbf{u}$ is the vorticity, $\nu$ is the kinematic viscosity and $\mathbf{U}_0$ is a mean flow velocity. $\mathbf{U}_s$ is the velocity of the solid at position $\mathbf{x}$ and at time $t$,

$$ \mathbf{U}_s(\mathbf{x}, t) = \dot{\mathbf{x}}_{cg}(t) + \dot{\alpha}_{cg}(t)[\mathbf{x} - \mathbf{x}_{cg}(t)]^{\perp}, \quad (31) $$

with $\mathbf{x}_{cg}$ being the center-of-gravity position vector of the solid, $\dot{\alpha}_{cg}$ being the angle of incidence. The dot denotes the derivative with respect to time, and $\mathbf{x}^{\perp} = (-y, x)$. The motion of the obstacle implies changing in time its mask function $\chi$. For its translation we rotate the phase of its Fourier coefficients:

$$ \chi(\mathbf{x}, t) = \chi_0(\mathbf{x} - \delta\mathbf{x}) \iff \tilde{\chi}(\mathbf{k}, t) = e^{-i\mathbf{k} \cdot \delta\mathbf{x}} \tilde{\chi}_0(\mathbf{k}), \quad (32) $$

where $\mathbf{k}$ is the wave number. Following [23], the solid body rotation at an angle $\beta$ is decomposed into three skewing operations:

$$ R(\beta) = \begin{bmatrix} \cos \beta & -\sin \beta \\ \sin \beta & \cos \beta \end{bmatrix} = \begin{bmatrix} 1 & -\tan(\beta/2) \\ 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 \\ \sin \beta & 1 \end{bmatrix} \begin{bmatrix} 1 & -\tan(\beta/2) \\ 0 & 1 \end{bmatrix}. \quad (33) $$

The mask function is smoothed to avoid Gibbs oscillations. This is done by solving the heat equation with a discontinuous mask function taken as initial condition. For details, we refer again to [23].

In order to include the time evolution of the passive scalar in such a configuration, the advection-diffusion equation is considered and is modified with a volume penalization method to impose zero scalar flux at the wall:

$$ \frac{\partial \theta}{\partial t} + (1 - \chi)(\mathbf{u} + \mathbf{U}_0)\nabla\theta + \chi(\mathbf{U}_s\nabla\theta) = \nabla \cdot ([\kappa(1 - \chi) + \eta_0\chi]\nabla\theta), \quad (34) $$

where $\theta$ is the passive scalar, $\kappa$ is the scalar diffusivity, and $\eta_0$ is the penalization parameter for the advection-diffusion of the passive scalar. Note the additional term in Eq. (34) compared to the original equation (Section 2): $\chi(\mathbf{U}_s\nabla\theta)$. This term is important because it transports the artificial scalar field, which is created inside the moving solid body to ensure the no-flux condition. Otherwise these numerical artifacts would be ejected inside the fluid domain once the solid left its position. As for the vorticity equation, we used a smoothed mask function to avoid Gibbs oscillations.

To solve the Eqs. (30) and (34), we use again a classical Fourier pseudo-spectral scheme with semi-implicit time integration, an Adams-Bashforth scheme for the nonlinear and the penalization terms and Euler backwards for the diffusion term in (30), which are both of second order. Note that in (34), the diffusion term is discretized explicitly with Adams-Bashforth.

### 4.1.2. Numerical test

In this part, we describe the test case used for validation. It consists in comparing a fixed obstacle in a mean flow $\mathbf{U}_0 = (0, 1)$ and a moving obstacle which moves in the fluid domain with a velocity equal to $\mathbf{U}_s = (0, -1)$. Thus we obtain two configurations:

* configuration 1: 'fixed' with $\mathbf{U}_0 = (0, 1)$ and $\mathbf{U}_s = (0, 0)$.

* configuration 2: 'moving' with $\mathbf{U}_0 = (0, 0)$ and $\mathbf{U}_s = (0, -1)$,

which are equivalent in the continuous formulation due to Galilean invariance.

For both configurations, the resolution of the computation is equal to $N = 512^2$ grid points and the time step is $\Delta t = 10^{-4}$. The penalization parameters are equal to $\eta_u = 10^{-4}$ and $\eta_0 = 10^{-7}$ for the Eqs. (30) and (34), respectively. The kinematic viscosity is set to $\nu = 10^{-2}$. The initial condition for the passive scalar is a circular blob placed in the domain at $(x, y) = (\pi, \pi/3)$. Its characteristic radius is $R = 0.7$. The Schmidt number $Sc = \nu/\kappa$ is unity and the domain size is $[0, 2\pi] \times [0, 2\pi]$.

Figs. 13 and 14 show visualizations of vorticity, $\omega = \nabla \times \mathbf{u}$, and the passive scalar field, $\theta$, for both configurations. Qualitatively, we observe the same results for the vorticity and the scalar. Moreover, for the scalar, oscillations with negatives values are present inside the solid which illustrates the numerical artifact in the solid due to the penalization method. As these values could be non-negligible, all the statistics of the scalar fluid are computed inside the fluid domain only and not in the full domain.
Fig. 15 presents the time evolution of the kinetic energy $E = 1/2\langle u^2 \rangle$ where $\langle \cdot \rangle$ is the space average and of the entrophy $Z = 1/2\langle \omega^2 \rangle$ for both configurations (proportional to the energy dissipation). The curves of the two configurations collapse for both, energy and entrophy, which gives us a first indication on the validity of the method used to compute the flow around a moving solid.

Fig. 16 shows the scalar variance $\sigma^2(t) = \langle(\theta(t) - \langle\theta(t)\rangle_x)^2\rangle_\gamma$, the scalar dissipation rate $\Lambda(t) = \langle|\nabla\theta(t)|^2\rangle_\gamma$, and the scalar gradient production $\Gamma(t) = \int_{\Omega_\gamma} \gamma(\mathbf{x}, t)d\mathbf{x}$ where the local scalar gradient production is $\gamma' = 1/2|\nabla\theta|\mathrm{scos}2\alpha$ and $\alpha$ is the angle between the scalar gradient and the eigenvector of the strain vector $s = \sqrt{s_1^2 + s_2^2}$ with $s_1 = \partial_x u_1 - \partial_y u_2$, see, e.g., [31]. Note that $\langle \cdot \rangle_\gamma$ denotes the spatial average computed only in the fluid domain. For both configurations, all three quantities do perfectly agree, which validates the method for moving solid obstacles, described by Eq. (34).

In the exact, continuous formulation the two cases are identical up to a change of the reference frame. However, the corresponding discretized problems are not algebraically equivalent. We show now that both approximated solutions converge to the same limit, which is taken from the fixed-cylinder case and which is assumed to be the reference solution of the continuous problem. The numerical simulations are performed with three different spatial resolutions: $N^2 = 256^2$, $512^2$ and $1024^2$.

Fig. 17 displays the difference between the 'fixed' and 'moving' cases for simulations with $N^2 = 1024^2$. At $t=0$ the fixed and the moving cylinder are not in the same position with respect to the grid. Therefore the vorticity field of the moving cylinder is interpolated on the grid of the fixed cylinder, and then the point-wise difference is calculated. This is shown in the left panel of Fig. 17. The right panel shows a similar quantity calculated for the passive scalar. The maximum error is localized near the solid for the vorticity and inside the solid for the scalar. Indeed for the scalar, the solution adapts itself inside the solid to insure Neumann conditions and these numerical artifacts are advected with the solid, which explains why the errors are concentrated inside the solid.

Fig. 18 displays the $L^\infty$ and $L^2$-norms of the difference of the vorticity fields for different resolutions. The values in the left panel are calculated in the whole computational domain, containing thus both, the solid and the fluid. In the right panel the contribution of the solid domain is excluded. The results show that the vorticity fields differ primarily inside the fluid domain. The difference is converging with the grid size with $N^{-1}$ for both the $L^2$ and $L^\infty$-norms.

In contrast, the errors for the passive scalar, shown in Fig. 19, are only converging inside the fluid domain. The lack of convergence inside the solid is explained by the different treatment of the advection term in the two cases. For the physical interpretation of the results this is insignificant as only the scalar inside the fluid domain is considered.

## 4.2. Results for a dynamical mixer

In this part we present numerical simulations of flow with passive scalar in a simplified dynamical mixing device. The study of such flows is relevant for many industrial applications for control and optimization, for example in chemical

<figure><img src="image_14.png" /><figcaption>Fig. 13. (Color online) Snapshots of vorticity at different time instants for configuration 'fixed' (top) and 'moving' (bottom).</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 14. Snapshots of the scalar fields at different time instants for the configuration 'fixed' (top) and 'moving' (bottom).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15. Time evolution of energy (left) and enstrophy (right) for both configurations.</figcaption></figure>

reactors. As penalization methods are very flexible concerning the choice of geometries, these techniques could have some impact on the optimization of industrial devices.

As computational domain, we consider a square with side length $2\pi$ and impose periodic boundary conditions. The fluid domain is confined inside a disk of radius $R = 2.8$. The cross shaped rotor of the mixer is located in the center of the domain and has four rectangular blades of size $L = 1$ and width $W = 0.05$. At $t = 0$ the fluid is at rest and the rotation of the blades starts impulsively with angular velocity of $\dot{\alpha}_{cg} = 1$. The initial condition of the passive scalar is a Gaussian blob with a characteristic size equal to $R_b = 0.79$. The blob is placed inside the domain at $(x,y) = (\pi - 1.5, \pi)$. The viscosity and the diffusivity are $\nu = \kappa = 10^{-4}$. The resolution is $N = 1024^2$.

Figs. 20 and 21 show respectively snapshots of vorticity and scalar fields for different time instants. The blades of the rotor create vortices which are ejected from the blades because of the detachment of the boundary layer. Negative vortices are created on the tips of the mixer and positive vortices are formed on the blades. Sometimes, two counter-rotative vortices form a dipole which advances and crashes into the wall. We can observe that the scalar is advected by the mean rotation and rolled-up by small vortices created by the blades which thus enhance the mixing.
<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Time evolution of scalar variance (top, left), scalar dissipation (top, right) and scalar production (bottom) for both configurations.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 17.** (Color online) Illustration of the pointwise error between the 'fixed' and 'moving' configurations. Difference between the vorticity fields (left) and the scalar fields (right) computed for the fixed and the moving cylinder at $t = 10$ for $N^2 = 1024^2$.</figcaption></figure>

Fig. 22 shows the time evolution of energy and entrophy. The energy increases rapidly with oscillations and then fluctuates around a constant value, while the entrophy oscillates even more and then decays non-monotonically due to the vorticity production at the boundaries.

Fig. 23 gives the time evolution of scalar variance and scalar dissipation rate. The scalar variance exhibits a strong decrease for small times while the decrease slows down for long times. The scalar dissipation rate confirms the previous figures, since this curve yields information on the slope of the scalar variance, it increases for small times and reaches a plateau, and then, for long times, decreases and reaches another plateau.
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 18. $L^\infty$ and $L^2$-norms of the error for the vorticity field between the 'fixed' and 'moving' configurations as a function of the number of grid points $N$ for the whole domain $\Omega$ (left) and the fluid domain $\Omega_f$ (right) at $t = 10$.</figcaption></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 19. $L^2$ and $L^\infty$-error for the scalar fields, computed in the whole domain (left) and computed in the fluid domain only (right) at $t = 10$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 20. (Color online) Snapshots of vorticity fields for the mixer.</figcaption></figure>

In conclusion, this example illustrates the flexibility and the potential of the volume penalization method for studying turbulent mixing in complex and time varying geometries.

## 5. Conclusions

The volume penalization method is an immersed boundary method which enables to impose efficiently in Fourier pseudo-spectral codes not only Dirichlet, but also Neumann boundary conditions. For the Neumann case we have shown
<figure><img src="image_24.png" /><figcaption>Fig. 21. Snapshots of scalar fields for the mixer.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 22. Time evolution of the energy (left) and enstrophy (right).</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 23. Time evolution of the scalar variance (left), the scalar dissipation (right).</figcaption></figure>

theoretically and numerically, considering simple examples, that the exact solution of the penalized equation converges with $\sqrt{\eta_0}$ to the exact solution of the equation satisfying Neumann boundary conditions. The convergence rate is thus the same as for the volume penalization method imposing Dirichlet conditions [15]. However, in contrast to the Dirichlet case the main advantage is that the penalization parameter $\eta_0$ can be chosen as small as desired without limiting the time step of the numerical method while still using explicit time discretization. Thus the penalization error can be minimized without imposing a stability limit to the numerical scheme. A drawback is the low regularity of the exact solution of the penalized equation which determines the convergence rate of the spectral method. The considered examples showed convergence of order $O(N^{-1/2})$ and $O(N^{-3/4})$ in one and two space dimensions, respectively.
Nevertheless the method is robust and allows to simulate flows imposing no-slip walls for the fluid and no-flux conditions for passive scalars in complexly shaped domains which even may vary in time. The volume penalization method is an efficient way to implement Neumann boundary conditions into a Fourier pseudo-spectral solver, which is well optimized using parallel fast Fourier transforms and widely used for turbulent flows. In addition no linear systems have to be solved and no additional numerical diffusion is added by the Fourier spectral discretization. The extension of the penalization method to three-dimensional flows is straightforward and thus parallel computations of flows in complex domains are feasible using available codes on massively parallel computers.

The further development of the volume penalization method for electrically conducting fluids is currently in progress. First results for magnetohydrodynamic turbulence in confined domains are promising and can be found in [32].

## Acknowledgments

We thank W.J.T. Bos for fruitful discussions.

## Appendix A. Summary of the proof for the convergence and error estimates theorem

In the following we present a short sketch of the proof of Theorem 2.1 which states convergence and error estimates of the penalized solution for the unsteady advection-diffusion equation. Further details can be found in the manuscript [30].

The domain $\Omega \subset \mathbb{R}^d$ (with $d=2$ or 3) is considered to be an open bounded set, generally chosen convex and polygonal; its boundary is denoted by $\Gamma = \partial\Omega$. The boundaries of the disjoint subdomains, i.e., the fluid and solid domains, $\Omega_f$ and $\Omega_s$, respectively, are Lipschitz continuous.

We consider the following unsteady advection-diffusion equation

$$ \partial_t \theta + \mathbf{u} \cdot \nabla \theta = \kappa \nabla^2 \theta \quad \text{for } \mathbf{x} \in \Omega_f, t > 0 \qquad (A.1) $$

completed with a homogeneous Neumann boundary condition, $\kappa \nabla \theta \cdot \mathbf{n} = 0$, for $\mathbf{x} \in \partial\Omega_f$ (where $\mathbf{n}$ denotes the outer normal vector on $\partial\Omega_f$) and a suitable initial condition, $\theta(t=0) = \theta_0$. The velocity $\mathbf{u}$ is supposed to be given and satisfies the incompressible Navier-Stokes equations. In particular we have $\nabla \cdot \mathbf{u} = 0$. The above problem is extended into the larger domain $\Omega$ using the penalization technique, i.e., we consider Eq. (2), with penalization parameter $\eta_0$.

The convergence proof and the derivation of the error estimates of the penalized problem can be sketched as follows: The skew symmetry of the advection term, together with a weak formulation of the penalized problem in $H^1$ allow to show the coercivity of the elliptic operator. Then it can be verified that the assumptions of J.-L. Lions' theorem for abstract parabolic problems, see e.g. [Chap. 3.4 in J.-L. Lions, E. Magenes, Problèmes aux limites non homogènes et applications, vol. 1, Dunod (Paris), 1968.], are satisfied which ensures the existence and uniqueness of the solution of the weak penalized problem. The error estimate of the homogeneous Neumann condition ($\sqrt{\eta_0}$) can be obtained from the energy estimate of the weak formulation using Gronwall's inequality. The error estimate ($\sqrt{\eta_0}$) in the solid domain $\Omega_s$ can then be directly deduced.

For the error estimate in the fluid domain $\Omega_f$ the situation is more complicated. Here the difference between the weak penalized problem and the weak limit problem ($\eta_0 \to 0$) is considered, which yields a weak equation for the error in the fluid domain. Using a particular test function and the skew symmetry of the advection term it can be shown that the latter indeed vanishes. Then applying successively different inequalities, i.e., the duality inequality, the trace inequality, mean Poincaré's and Cauchy-Schwarz inequalities, the time evolution of the error can be estimated and then bounded by applying Young's inequality. Integration of the estimate in time and using that the error of the initial condition vanishes, yields, by using the error bound of the Neumann condition and by applying Gronwall's inequality, the desired error estimate of order $\sqrt{\eta_0}$ in $\Omega_f$.

## References

[1] J. Ferziger, M. Peric, Numerical Methods in Fluid Dynamics, Springer-Verlag, Berlin, 1996.

[2] J.F. Thompson, Z.U.A. Warsi, C.W. Mastin, Numerical Grid Generation: Foundations and Applications, Elsevier North-Holland, Inc, New York, NY, 1985.

[3] S.P. Spekreijse, Elliptic grid generation based on laplace equations and algebraic transformations, J. Comput. Phys. 118 (1995) 38-61.

[4] R. Courant, Variational methods for the solution of problems of equilibrium and vibrations, Bull. Amer. Math. Soc. 49 (1943) 1-23.

[5] R. Glowinski, Y. Kuznetsov, On the solution of the Dirichlet problem for linear elliptic operators by a distributed Lagrange multiplier method, C. R. Acad. Sci. Paris Sér. I 327 (1999) 693-698.

[6] Y. Cheny, O. Botella, The LS-STAG method: a new immersed boundary/level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties, J. Comput. Phys. 229 (2010) 1043-1076.

[7] C.S. Peskin, Flow patterns around heart valves: a digital computer method for solving the equations of motion. Ph.D. thesis, Albert Einstein Coll. Med., 1972.

[8] P. Angot, C.-H. Bruneau, P. Fabrie, A penalisation method to take into account obstacles in viscous flows, Numer. Math. 81 (1999) 497.

[9] C.S. Peskin, The immersed boundary method, Acta Numer. 11 (2002) 479-517.

[10] R. Mittal, G. Iaccarino, Immersed boundary methods, Annu. Rev. Fluid Mech. 37 (2005) 23961.

[11] P. Joly, L. Rhaouti, Fictitious domains, H(div) finite elements and Neumann conditions: the inf-sup condition, C. R. Acad. Sci. Paris Sér. I 328 (1999) 1225-1230.

[12] R. Glowinski, T.W. Pan, Error estimates for fictitious domain/penalty/finite elements methods, Calcolo 29 (1-2) (1992) 125-141.

[13] E. Arquis, J.-P. Caltagirone, Sur les conditions hydrodynamiques au voisinage d'une interface milieu fluide milieu poreux: application à la convection naturelle, C. R. Acad. Sci. Paris Sér. II 299 (1984) 14.

[14] P. Angot, A unified fictitious domain model for general embedded boundary conditions, C. R. Acad. Sci. Paris Sér. I 341 (2005) 683-688.

[15] G. Carbou, P. Fabrie, Boundary layer for a penalization method for viscous incompressible flow, Adv. Differ. Equat. 8 (2003) 1453-2480.
[16] N.K.-R. Kevlahan, J.-M. Ghidaglia, Computation of turbulent flow past an array of cylinders using a spectral method with Brinkman penalization, Eur. J. Mech. B/Fluids 20 (2001) 333-350.

[17] K. Schneider, Numerical simulation of the transient flow behaviour in chemical reactors using a penalisation method, Comput. Fluids 34 (2005) 1223-1238.

[18] H.J.H. Clercx, C.-H. Bruneau, The normal and oblique collision of a dipole with a no-slip boundary, Comput. Fluids 35 (2006) 245.

[19] G.H. Keetels, U. d'Ortona, W. Kramer, H.J.H. Clercx, K. Schneider, G.J.F. van Heijst, Fourier spectral and wavelet solvers for the incompressible Navier-Stokes equations with volume penalization: convergence of a dipole-wall collision, J. Comput. Phys. 227 (2007) 919-945.

[20] K. Schneider, M. Farge, Decaying two-dimensional turbulence in a circular container, Phys. Rev. Lett. 95 (2005) 244-502.

[21] B. Kadoch, W.J.T. Bos, K. Schneider, Extreme Lagrangian acceleration in confined turbulent flow, Phys. Rev. Lett. 100 (2008) 184-503.

[22] I. Ramière, P. Angot, M. Belliard, A fictitious domain approach with spread interface for elliptic problems with general boundary conditions, Comput. Methods Appl. Mech. Eng. 196 (2007) 766-781.

[23] D. Kolomenskiy, K. Schneider, A Fourier spectral method for the Navier-Stokes equations with volume penalization for moving solid obstacles, J. Comput. Phys. 228 (2009) 5687-5709.

[24] D.M. Hobbs, F.J. Muzzio, Reynolds number effects on laminar mixing in the Kenics static mixer, Chem. Eng. J. 70 (1998) 93-103.

[25] J.M. Ottino, The Kinematics of Mixing: Stretching, Chaos, and Transport, Cambridge University Press, Cambridge, 1989.

[26] Z. Warhaft, Passive scalars in turbulent flows, Ann. Rev. Fluid Mech. 32 (2000) 203-240.

[27] G. Boffetta, F. De Lillo, A. Mazzino, Peripheral mixing of passive scalar at small Reynolds number, J. Fluid Mech. 624 (2009) 151-158.

[28] F. Schwertfing, M. Manhart, DNS of passive scalar transport in turbulent channel flow at high Schmidt numbers, Int. J. Heat Fluid Flow 28 (2007) 1204-1214.

[29] A. Paccou, G. Chiavassa, J. Liandrat, K. Schneider, A penalization method applied to the wave equation, C. R. Méc. 333 (1) (2005) 79-85.

[30] P. Angot, B. Kadoch, D. Kolomenskiy, K. Schneider, Convergence analysis of a penalty method for scalar transport and mixing in Navier-Stokes flows, Manuscript 2012.

[31] C. Beta, K. Schneider, M. Farge, Wavelet filtering to study mixing in two-dimensional isotropic turbulence, Commun. Nonlin. Sci. Num. Sim. 8 (2003) 537-545.

[32] K. Schneider, S. Neffaa, W.J.T. Bos, A pseudo-spectral method with volume penalisation for magnetohydrodynamic turbulence in confined domains, Comput. Phys. Commun. 182 (1) (2011) 2-7.