1. Introduction
Immersed boundary techniques, including penalization approaches, are nowadays commonly employed to solve boundary or initial boundary value problems in complex geometries. They consist of embedding the original, possibly complex, spatial domain inside a bigger domain having a simpler geometry, for example a Cartesian geometry, while keeping the boundary conditions approximately enforced thanks to new terms that are added to the equations. Historically, the penalization technique can be traced back to Courant (Reference Courant1943) in the context of constrained optimization. Later, Saul’ev (Reference Saul’ev1963) applied it in the context of fictitious domain methods for immersed boundaries. In a classical paper (Peskin Reference Peskin1977) an immersed boundary method was proposed for computing the blood flow in a beating heart. Some history on immersed boundary techniques can be found in Angot et al. (Reference Angot, Kadoch, Kolomenskiy and Schneider2014), and for detailed reviews we refer to the classical review papers Peskin (Reference Peskin2002), Mittal & Iaccarino (Reference Mittal and Iaccarino2005).
In the current work we focus on one particular example, the volume penalization method (Angot, Bruneau & Fabrie Reference Angot, Bruneau and Fabrie1999) which, inspired by the physical intuition that a solid wall is similar to a vanishingly porous medium, uses the Brinkman–Darcy drag force as penalization term. The main advantage of such penalized equations is that they can be discretized independently of the geometry of the original problem, since the latter has been encoded into the penalization terms. Such a simplification permits a massive reduction in solver development time, since it avoids the issues associated with the design and management of the grid, allowing for example the use of simple fast Fourier transform (FFT)-based spectral solvers in Cartesian geometries. Using spectral solvers no linear systems have to be solved, e.g., to impose the incompressibility of the magnetic and velocity fields. Furthermore, the penalized equations are solved without introducing additional discretization errors, and thus effects of numerical diffusion and dispersion present in low-order numerical schemes are avoided. Moreover, parallel FFT libraries (e.g. P3DFFT) are available on the current supercomputers. The advantage of the penalization method becomes even more substantial when the geometry is time-dependent, as in the case of moving obstacles, or when fluid–structure interaction is taken into account.
In the following we give a non-exhaustive overview of the development for different applications using volume penalization. Two-dimensional hydrodynamic flows in staggered or inline tube bundles have been computed in Schneider (Reference Schneider2005) and a dipole–wall collision benchmark using different numerical discretizations of the penalized Navier–Stokes equations has been obtained in Keetels et al. (Reference Keetels, d’Ortona, Kramer, Clercx, Schneider and van Heijst2007). Final states in two-dimensional decaying hydrodynamic turbulence in different geometries have been studied in Schneider & Farge (Reference Schneider and Farge2008) and flows past flat plates focusing on the influence of the plate’s geometry in Schneider et al. (Reference Schneider, Paget-Goy, Verga and Farge2014). Nguyen van yen, Farge & Schneider (Reference Nguyen van yen, Farge and Schneider2011) computed dipole–wall collisions and derived the scaling of the energy dissipation in the large Reynolds number limit. They showed that it converges indeed to a finite value. The penalized wave equation has been analysed in Paccou et al. (Reference Paccou, Chiavassa, Liandrat and Schneider2005). Applications to moving obstacles and fluid–structure interaction problems can be found in Kolomenskiy & Schneider (Reference Kolomenskiy and Schneider2009), Kolomenskiy et al. (Reference Kolomenskiy, Moffatt, Farge and Schneider2011), Kolomenskiy, Engels & Schneider (Reference Kolomenskiy, Engels and Schneider2013) and Engels et al. (Reference Engels, Kolomenskiy, Schneider and Sesterhenn2013). In Kreuzahler et al. (Reference Kreuzahler, Schulz, Homann, Ponty and Grauer2014) a numerical study using the volume penalization method for impeller-driven von Kármán flows has been performed. Different simulations of compressible flows in complex geometries using volume penalization have been described in Liu & Vasilyev (Reference Liu and Vasilyev2007) and Boiron, Chiavassa & Donat (Reference Boiron, Chiavassa and Donat2009). An adaptive penalty method for incompressible Navier–Stokes has been proposed by Shirokoff & Nave (Reference Shirokoff and Nave2015).
The extension of penalization techniques to magnetohydrodynamic (MHD) flows is more recent. Simulations of the self-organization of confined plasma in two space dimensions studying the influence of the confinement geometry have been presented in Bos, Neffaa & Schneider (Reference Bos, Neffaa and Schneider2008) and Neffaa, Bos & Schneider (Reference Neffaa, Bos and Schneider2008). The numerical method has been extended to three space dimensions and is described and benchmarked in detail in Morales et al. (Reference Morales, Leroy, Bos and Schneider2014b ). The spontaneous spin-up of plasma in toroidal geometries using the viscous-resistive MHD equations has been discovered in Morales et al. (Reference Morales, Bos, Schneider and Montgomery2012). The effect of toroidicity in reversed field pinch (RFP) devices is studied in Morales et al. (Reference Morales, Bos, Schneider and Montgomery2014a ). Roberts et al. (Reference Roberts, Leroy, Morales, Bos and Schneider2014) investigated helically forced MHD flows in cylindrical geometries and showed the persistence of helical modes even when the dynamics becomes turbulent. A penalty method for hyperbolic partial differential equations modelling the edge plasma transport in a tokamak has been proposed in Angot, Auphan & Gues (Reference Angot, Auphan and Gues2014).
The aim of this paper is to review the volume penalization method and to illustrate its potential for applications in fluid and plasma turbulence in complex domains. Typically, Dirichlet boundary conditions are considered corresponding to imposing the value of the solution at the boundary. Extensions for dealing with Neumann boundary conditions, i.e. imposing values of the derivative of the solution at the boundary, have been proposed in Kadoch et al. (Reference Kadoch, Kolomenskiy, Angot and Schneider2012) and Kolomenskiy, Nguyen van yen & Schneider (Reference Kolomenskiy, Nguyen van yen and Schneider2015) based on previous work in Ramière, Angot & Belliard (Reference Ramière, Angot and Belliard2007). An alternative approach for imposing either Dirichlet or Neumann boundary conditions which is based on sharp interface methods has been introduced in Mittal et al. (Reference Mittal, Dong, Bozkurttas, Najjar, Vargas and von Loebbecke2008). Considering simple examples in one space dimension allows an understanding and analysis of the convergence behaviour of the penalization techniques when the penalization parameter tends to zero. We show that, for a given numerical discretization of the penalized equations, there exists a value of the penalization parameter, corresponding to a balance between penalization and discretization errors, below which no further gain in precision is achieved. These results shed light on the behaviour of volume penalization schemes when solving the penalized equations, outline the limitations of the method, and give indications on how to choose the penalization parameter, at least in simple test cases. Nevertheless, for practical applications a series of computations may be necessary to determine the actual value of the parameters, as discussed e.g., in Engels et al. (Reference Engels, Kolomenskiy, Schneider and Sesterhenn2015). Finally, different illustrations will be given for hydro- and magnetohydrodynamic problems in the turbulent regime, including a study of the spatio-temporal self-organization of visco-resistive magnetohydrodynamics in toroidal geometries with different poloidal cross-sections while imposing curl-free toroidal magnetic and electric fields.
The outline of the paper is the following. First, we present a short primer on penalization. Then we describe in some detail the volume penalization to impose either Dirichlet or Neumann boundary conditions and we present some simple one-dimensional examples. Applications to fluid and plasma turbulence in two and three dimensions illustrate the properties of the method. Finally, we set out the conclusions together with some perspectives for future work.
2. A short primer on penalization
To illustrate the idea of volume penalization we consider a boundary (BVP) or initial boundary value problem (IBVP) for a partial differential equation (PDE) in a domain
${\it\Omega}_{f}\subset \mathbb{R}^{d}$
, written in the following abstract operator equation:
$$\begin{eqnarray}Lu=f\quad \text{for }\boldsymbol{x}\in {\it\Omega}_{f},\end{eqnarray}$$
completed with boundary conditions
$Bu=g$
at
$\partial {\it\Omega}_{f}$
and additional initial conditions in the case of an IBVP. The differential operator
$L$
stands, for example, for the Laplace operator,
$L=-{\rm\nabla}^{2}$
, or the Navier–Stokes or Maxwell operator. The boundary conditions can be of Dirichlet or Neumann type, e.g.,
$u=g$
or
$\partial u/\partial n=g$
, respectively, where
$n$
is the outer normal of the domain. Solving (2.1) numerically requires domain-fitted grids using, e.g. finite elements. The grid generation and the numerical solution of the resulting discretized problem can be demanding, see for example Ferziger & Peric (Reference Ferziger and Peric1996).
An alternative are penalization methods which embed the problem posed in the complexly shaped domain
${\it\Omega}_{f}$
into a larger simple domain (see figure 1), typically of rectangular shape
${\it\Omega}$
, i.e.,
${\it\Omega}_{f}\subset {\it\Omega}$
. The advantage is that fast solvers are available for such problems, which can furthermore be easily parallelized.
Sketch of the fluid domain
${\it\Omega}_{f}$
immersed in the solid domain
${\it\Omega}_{s}$
. The boundary of the fluid domain is denoted by
${\it\Sigma}=\partial {\it\Omega}_{f}$
. The total computational domain is
${\it\Omega}={\it\Omega}_{f}\cup {\it\Omega}_{s}$
.
In the Dirichlet case the penalized problem thus reads
$$\begin{eqnarray}Lu_{{\it\eta}}=f-\frac{1}{{\it\eta}}{\it\chi}(Bu-g)\quad \text{for }\boldsymbol{x}\in {\it\Omega}={\it\Omega}_{f}\cup {\it\Omega}_{s},\end{eqnarray}$$
where the boundary conditions have been included into the equation and an additional parameter, the small penalization parameter
${\it\eta}>0$
, has been introduced. All information about the geometry of
${\it\Omega}_{f}$
has been encoded into the mask function
${\it\chi}$
, which is defined as
$$\begin{eqnarray}{\it\chi}(\boldsymbol{x})=\left\{\begin{array}{@{}ll@{}}0\quad & \text{for }\boldsymbol{x}\in {\it\Omega}_{f}\\ 1\quad & \text{elsewhere}.\end{array}\right.\end{eqnarray}$$
In the domain
${\it\Omega}_{f}$
the mask function vanishes identically and thus the original PDE (2.1) is satisfied.
The difference between the solution
$u$
of the original PDE (2.1) and the solution of the penalized problem
$u_{{\it\eta}}$
(2.2) is the penalization (or modelling) error, which depends on the size of the penalization parameter
${\it\eta}$
,
$$\begin{eqnarray}\Vert u-u_{{\it\eta}}\Vert \propto {\it\eta}^{{\it\alpha}},\end{eqnarray}$$
where
${\it\alpha}$
describes the order of the penalization method and
$\Vert \cdot \Vert$
is a suitable norm, e.g. the energy norm
$\Vert f\Vert _{2}=(\int _{{\it\Omega}_{f}}|f(\boldsymbol{x})|^{2}\,\text{d}\boldsymbol{x})^{1/2}$
or the
$L^{1}$
or
$L^{\infty }$
norm. For convergence of the method the solution of the penalized problem
$u_{{\it\eta}}$
should tend towards the solution of the original problem
$u$
, i.e.
$\lim _{{\it\eta}\rightarrow 0}\Vert u-u_{{\it\eta}}\Vert \rightarrow 0$
. Thus
${\it\alpha}>0$
is required for convergence of the method. For the volume penalization we have
${\it\alpha}=1/2$
(Angot et al.
Reference Angot, Bruneau and Fabrie1999; Carbou & Fabrie Reference Carbou and Fabrie2003).
Applying a numerical method to the penalized PDE (2.2) we obtain the discretized penalized equation
$$\begin{eqnarray}L^{N}u_{{\it\eta}}^{N}=f^{N}-\frac{1}{{\it\eta}}{\it\chi}^{N}(Bu-g)\quad \text{for }\boldsymbol{x}\in {\it\Omega},\end{eqnarray}$$
with
${\rm\Delta}x\propto 1/N$
, where
$N$
denotes the number of grid points and
$L^{N}$
is the discretized version of the operator
$L$
. In the case of an IBVP a suitable time discretization is also necessary. Explicit time discretization of the penalization term implies a time step restriction for stability. Typically, the time step is limited by
${\it\eta}$
and hence the problem becomes stiff in the limit
${\it\eta}\rightarrow 0$
.
The solution of the discretized penalized equations
$u_{{\it\eta}}^{N}$
, given that the discretization is consistent and the numerical scheme is stable, converges towards the exact solution
$u_{{\it\eta}}$
of the penalized equation (2.2). The discretization error,
$$\begin{eqnarray}\Vert u_{{\it\eta}}-u_{{\it\eta}}^{N}\Vert \propto \left(\frac{1}{N}\right)^{{\it\beta}},\end{eqnarray}$$
depends not only on the order of the numerical scheme, but is also limited by the regularity, i.e. the smoothness of the exact solution
$u_{{\it\eta}}$
of the penalized problem. The order
${\it\beta}$
is thus determined as the mimum of the order of the underlying numerical scheme and the regularity of the exact penalized solution.
Finally, the question in computations using the penalization method concerns the error between the exact solution
$u$
of the BVP (or IBVP) given by (2.1) and the numerical solution of the penalized equation (2.5). This total error has two contributions: the modelling error which is due to the penalization; and the discretization error which is due to the numerical solution of the discretized problem. Applying the triangle inequality the following estimate holds:
$$\begin{eqnarray}\Vert u-u_{{\it\eta}}^{N}\Vert \leqslant \Vert u-u_{{\it\eta}}\Vert +\Vert u_{{\it\eta}}-u_{{\it\eta}}^{N}\Vert ,\end{eqnarray}$$
and thus we obtain a bound for the total error. A straightforward argument would suggest choosing very small values for
${\it\eta}$
to minimize the modelling error. However this would imply that the problem becomes stiff, one loses regularity in the exact solution of the penalized problem and the leading-order constant of the discretization error blows up and a very fine grid would become necessary. Thus to optimize the estimate both errors should be of the same order of magnitude. This shows that the penalization parameter
${\it\eta}$
and the numerical resolution
$N$
of the discretized problem are coupled and should be chosen accordingly. For further discussion we refer to Nguyen van yen, Kolomenskiy & Schneider (Reference Nguyen van yen, Kolomenskiy and Schneider2014).
3. Volume penalization for Dirichlet and Neumann boundary conditions
For illustration we consider first a simple toy problem, the Poisson equation in one space dimension,
$$\begin{eqnarray}-\frac{\text{d}^{2}u}{\text{d}x^{2}}=f\quad \text{for }x\in {\it\Omega}_{f},\end{eqnarray}$$
completed with either homogeneous Dirichlet (
$u=0$
for
$x\in \partial {\it\Omega}_{f}$
) or Neumann boundary conditions (
$\text{d}u/\text{d}x=0$
for
$x\in \partial {\it\Omega}_{f}$
). We choose
${\it\Omega}_{f}= ]0,{\rm\pi}[$
and for the right-hand side
$f$
a trigonometric function, i.e.,
$f(x)=m^{2}\sin mx$
and
$f(x)=m^{2}\cos mx$
, respectively, with
$m\in \mathbb{N}$
. The exact solution is
$u(x)=\sin mx$
and
$u(x)=\cos mx+C$
in the Dirichlet and Neumann case, respectively. The additive constant
$C\in \mathbb{R}$
shows that the solution is not unique for Neumann boundary conditions. To guarantee in this case the existence of a solution, the right-hand side
$f$
has to satisfy the compatibility condition,
$\int _{0}^{{\rm\pi}}f(x)\,\text{d}x=u^{\prime }(x={\rm\pi})-u^{\prime }(x=0)=0$
.
3.1. The penalized Poisson equation
Now we replace the original problems by the penalized problems. The domain
${\it\Omega}_{f}= ]0,{\rm\pi}[$
is embedded in the larger domain
${\it\Omega}= ]0,2{\rm\pi}[$
. Thus we have
${\it\Omega}={\it\Omega}_{f}\cup {\it\Omega}_{s}$
, where
${\it\Omega}_{s}$
is the penalization domain. For further simplification we impose periodic boundary conditions at the boundary of
${\it\Omega}$
.
In the Dirichlet case we obtain the penalized Poisson equation
$$\begin{eqnarray}-\frac{\text{d}^{2}u_{{\it\eta}}}{\text{d}x^{2}}+\frac{1}{{\it\eta}}\;{\it\chi}\,u_{{\it\eta}}=f\quad \text{for }x\in ]0,2{\rm\pi}[,\end{eqnarray}$$
while in the Neumann case we have
$$\begin{eqnarray}-\frac{\text{d}}{\text{d}x}((1-{\it\chi})+{\it\eta}{\it\chi})\frac{\text{d}}{\text{d}x}u_{{\it\eta}}=f\quad \text{for }x\in ]0,2{\rm\pi}[.\end{eqnarray}$$
The penalization term in the Dirichlet case forces the solution to vanish inside the penalization domain and thus imposes homogeneous Dirichlet conditions at the interface. In the Neumann case the boundary condition corresponds to zero flux through the interface of
${\it\Omega}_{f}$
and
${\it\Omega}_{s}$
. This can be achieved by imposing vanishing diffusivity (of order
${\it\eta}$
) inside the penalization domain. The function
$f$
is extended in the larger domain by zero padding, i.e.,
$f(x)=0$
for
$x\in {\it\Omega}_{s}$
. Again both problems can be solved analytically in each sub-domain and we obtain in the Dirichlet case (Nguyen van yen et al.
Reference Nguyen van yen, Kolomenskiy and Schneider2014)
$$\begin{eqnarray}u_{{\it\eta}}(x)=\left\{\begin{array}{@{}l@{}}\sin mx+A_{1}x+A_{2}\quad \text{for }x\in ]0,{\rm\pi}[\quad \\ \displaystyle \frac{m^{2}{\it\eta}}{1+{\it\eta}m^{2}}\sin mx+B_{1}\exp (-x/\sqrt{{\it\eta}})+B_{2}\exp (x/\sqrt{{\it\eta}})\quad \text{for }x\in ]{\rm\pi},2{\rm\pi}[,\quad \end{array}\right.\end{eqnarray}$$
and in the Neumann case (Kolomenskiy et al. Reference Kolomenskiy, Nguyen van yen and Schneider2015)
$$\begin{eqnarray}u_{{\it\eta}}(x)=\left\{\begin{array}{@{}ll@{}}\cos mx+A_{1}x+A_{2}\quad & \text{for }x\in ]0,{\rm\pi}[\\ B_{1}x+B_{2}\quad & \text{for }x\in ]{\rm\pi},2{\rm\pi}[.\end{array}\right.\end{eqnarray}$$
Imposing continuity of the solution and of the derivative at
$x=0=2{\rm\pi}$
and
$x={\rm\pi}$
the coefficients (
$A_{1},A_{2},B_{1},B_{2}$
) can be determined. The corresponding values can be found in Kolomenskiy et al. (Reference Kolomenskiy, Nguyen van yen and Schneider2015) and Nguyen van yen et al. (Reference Nguyen van yen, Kolomenskiy and Schneider2014). Note that in the Neumann case only three out of the four coefficients can be determined as the solution is not unique. The penalization error
$\Vert u-u_{{\it\eta}}\Vert$
can thus be explicitly computed and we find in the Dirichlet case the expected
$O(\sqrt{{\it\eta}})$
behaviour, while in the Neumann case an
$O({\it\eta})$
behaviour is found which is better than the
$O(\sqrt{{\it\eta}})$
shown in Kadoch et al. (Reference Kadoch, Kolomenskiy, Angot and Schneider2012) for the heat equation. In figure 2(a,b) the exact solution of the penalized Dirichlet problem is plotted for two values of
${\it\eta}$
(figure 2
a) together with the first derivative (figure 2
b). We find that for decreasing
${\it\eta}$
the penalized solution tends towards the solution of the unpenalized problem. We also observe the existence of a boundary layer in the penalized domain (at the interface of
${\it\Omega}_{f}$
and
${\it\Omega}_{s}$
, close to
$x={\rm\pi}$
) which becomes steeper when
${\it\eta}$
becomes smaller. For
${\it\eta}=10^{-6}$
we find that the first derivative becomes almost discontinuous at
$x={\rm\pi}$
. This shows that the regularity of
$u_{{\it\eta}}$
is lost in the limit of small
${\it\eta}$
and the problem becomes stiff.
(a) Exact solution
$u_{{\it\eta}}$
and (b) its first derivative
$u_{{\it\eta}}^{\prime }$
of the penalized Poisson equation with Dirichlet boundary conditions for
${\it\eta}=10^{-2}$
and 10−6. Convergence of the second-order finite difference scheme. (c) Error with respect to the exact solution of the Dirichlet problem in
${\it\Omega}_{f}$
. From Nguyen van yen et al. (Reference Nguyen van yen, Kolomenskiy and Schneider2014). Reproduced with kind permission of Springer Science and Business Media.
3.2. The discretized penalized Poisson equation
To solve the penalized equations numerically we apply a second-order finite difference discretization to the penalized problems (3.2) and (3.3). The computational domain
${\it\Omega}=[0,2{\rm\pi}]$
is discretized with
$N$
grid points
$x_{i}=i/(2{\rm\pi}),~i=0,\ldots ,N-1$
, applying periodic boundary conditions. We obtain the following linear system in the Dirichlet case:
$$\begin{eqnarray}\left(-\unicode[STIX]{x1D63F}_{2}+\frac{1}{{\it\eta}}{\bf\chi}\unicode[STIX]{x1D644}\right)\boldsymbol{U}=\boldsymbol{F},\end{eqnarray}$$
with the vectors
$\boldsymbol{U}=(u(x_{0}),\ldots ,u(x_{N-1}))$
,
${\bf\chi}=({\it\chi}(x_{0}),\ldots ,{\it\chi}(x_{N-1}))$
,
$\boldsymbol{F}=(f(x_{0}),\ldots ,f(x_{N-1}))$
in
$\mathbb{R}^{N}$
, the identity matrix
$\unicode[STIX]{x1D644}$
, and where
$$\begin{eqnarray}\unicode[STIX]{x1D63F}_{2}=\frac{1}{h^{2}}\left(\begin{array}{@{}ccccc@{}}-2 & 1 & & & 1\\ 1 & -2 & 1 & & \\ & & & \ddots & \\ & & 1 & -2 & 1\\ 1 & & & 1 & -2\end{array}\right)\end{eqnarray}$$
is the second-order central finite difference operator with
$h=2{\rm\pi}/N$
. The resulting tridiagonal system can be solved using standard numerical linear algebra tools.
In the Neumann case we get a second-order approximation with the following finite difference scheme:
$$\begin{eqnarray}-{\textstyle \frac{1}{2}}(\unicode[STIX]{x1D63F}_{F}{\it\Theta}\unicode[STIX]{x1D63F}_{B}+\unicode[STIX]{x1D63F}_{B}{\it\Theta}\unicode[STIX]{x1D63F}_{F})\,\boldsymbol{U}=\boldsymbol{F},\end{eqnarray}$$
where
${\it\Theta}=[{\it\theta}(x_{0}),{\it\theta}(x_{1}),\ldots ,{\it\theta}(x_{N-1})]$
with
${\it\theta}(x_{i})=(1-{\it\chi}(x_{i})+{\it\eta}{\it\chi}(x_{i}))$
. The first-order backward and forward finite difference operators
$\unicode[STIX]{x1D63F}_{B}$
and
$\unicode[STIX]{x1D63F}_{F}$
are defined, respectively, by
$$\begin{eqnarray}\unicode[STIX]{x1D63F}_{B}=\frac{1}{h}\left(\begin{array}{@{}ccccc@{}}1 & & & & -1\\ -1 & 1 & & & \\ & & & \ddots & \\ & & & -1 & 1\end{array}\right),\quad \unicode[STIX]{x1D63F}_{F}=\frac{1}{h}\left(\begin{array}{@{}ccccc@{}}-1 & 1 & & & \\ & -1 & 1 & & \\ & & & \ddots & \\ 1 & & & & -1\end{array}\right).\end{eqnarray}$$
The linear system (3.8) is singular, the matrix has an eigenvalue 0 and a solution only exists if the right-hand side
$\boldsymbol{F}$
is in its image. Special care has to be taken in solving the linear system using either the pseudo-inverse, or removing one equation.
Performing numerical experiments Kolomenskiy et al. (Reference Kolomenskiy, Nguyen van yen and Schneider2015), Nguyen van yen et al. (Reference Nguyen van yen, Kolomenskiy and Schneider2014) showed that second-order convergence of the numerical solution of the penalized equations towards the exact solution can be obtained, under the condition that
${\it\eta}$
is sufficiently small. Moreover, for fixed
${\it\eta}$
the total error (modelling plus discretization error) has a minimum at
$N\propto 1/\sqrt{{\it\eta}}$
. Figure 2(c) illustrates this for the Dirichlet case. Varying
${\it\eta}$
with
$1/N$
like
${\it\eta}\propto 1/N^{2}$
thus yields the best results in terms of minimum total error (Nguyen van yen et al.
Reference Nguyen van yen, Kolomenskiy and Schneider2014). Nevertheless it must be mentioned that for spectral or fourth-order finite difference discretizations only first-order convergence can be observed. The explanation is that the linear system using second-order finite differences is exactly the same for the penalized problem as for the original Dirichlet problem at all points inside the fluid domain, except at
$x_{1}$
and
$x_{N/2-1}$
where values of
$O(\sqrt{{\it\eta}})$
are involved (Nguyen van yen et al.
Reference Nguyen van yen, Kolomenskiy and Schneider2014). Hence for sufficiently small values of
${\it\eta}$
second-order convergence can be found. Note that for fourth-order finite differences the stencil is larger, thus more grid points are affected by the penalization term and the second-order convergence is lost. Spectral discretizations of the penalized operators have been analyzed in Nguyen van yen et al. (Reference Nguyen van yen, Kolomenskiy and Schneider2014) motivated by the work of Min & Gottlieb (Reference Min and Gottlieb2003), where Fourier approximations of elliptic problems with discontinuous coefficients have been studied.
3.3. The penalized Navier–Stokes and Maxwell equations
Analogously to the penalized Poisson equation we can impose no slip and no penetration, i.e., Dirichlet boundary conditions, in the incompressible Navier–Stokes equations by adding a penalization term to the momentum equation,
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial \boldsymbol{u}}{\partial t}+{\bf\omega}\times \boldsymbol{u}+\boldsymbol{{\rm\nabla}}{\it\Pi}-{\it\nu}{\rm\nabla}^{2}\boldsymbol{u}=-\frac{1}{{\it\eta}}{\it\chi}(\boldsymbol{u}-\boldsymbol{u}_{wall}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$
${\bf\omega}=\boldsymbol{{\rm\nabla}}\times \boldsymbol{u}$
is the vorticity and
${\it\Pi}=p+\boldsymbol{u}^{2}/2$
is the modified pressure which satisfies
$-{\rm\nabla}^{2}{\it\Pi}=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\bf\omega}\times \boldsymbol{u}+(1/{\it\eta}){\it\chi}(\boldsymbol{u}-\boldsymbol{u}_{wall}))$
. The kinematic viscosity is denoted by
${\it\nu}$
. The velocity at the wall,
$\boldsymbol{u}_{wall}$
, does vanish identically in the case of fixed walls, or is prescribed in the case of moving walls. In the limit
${\it\eta}\rightarrow 0$
the solution of the penalized Navier–Stokes equations tends to the solution of the Navier–Stokes equations with no-slip boundary conditions and the penalization error is of
$O(\sqrt{{\it\eta}})$
, as shown in Angot et al. (Reference Angot, Bruneau and Fabrie1999) and Carbou & Fabrie (Reference Carbou and Fabrie2003).
The induction equation which describes the evolution of the magnetic field
$\boldsymbol{B}$
can be penalized in a similar way and we obtain
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial \boldsymbol{B}}{\partial t}-\boldsymbol{{\rm\nabla}}\times (\boldsymbol{u}\times \boldsymbol{B})-{\it\lambda}{\rm\nabla}^{2}\boldsymbol{B}=-\frac{1}{{\it\eta}}{\it\chi}(\boldsymbol{B}-\boldsymbol{B}_{wall}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{B}=0, & \displaystyle\end{eqnarray}$$
${\it\lambda}$
is the magnetic diffusivity. The magnetic field at the wall,
$\boldsymbol{B}_{wall}$
, can be freely chosen and hence not all components have to be penalized. For example choosing
$\boldsymbol{B}_{wall}=B_{\Vert }$
, where
$B_{\Vert }$
is the component of
$\boldsymbol{B}$
parallel to the wall, only penalizes the normal component and leaves the parallel component free. This allows modelling perfectly conducting boundary conditions. Note that in the MHD case no mathematically rigorous convergence theorem has been proven yet, but in Morales et al. (Reference Morales, Leroy, Bos and Schneider2014b
) asymptotic arguments for estimating the penalization error have been given.
For the transport of a passive scalar
${\it\xi}$
we consider the advection–diffusion equation where no-flux conditions
$\boldsymbol{{\rm\nabla}}{\it\xi}\boldsymbol{\cdot }n=0$
, i.e., homogeneous Neumann boundary conditions, are imposed (Kadoch et al.
Reference Kadoch, Kolomenskiy, Angot and Schneider2012),
$$\begin{eqnarray}\frac{\partial {\it\xi}}{\partial t}+[(1-{\it\chi})\boldsymbol{u}]\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\xi}=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }[{\it\kappa}(1-{\it\chi})+{\it\eta}_{{\it\xi}}{\it\chi}]\boldsymbol{{\rm\nabla}}{\it\xi},\end{eqnarray}$$
where
${\it\eta}_{{\it\xi}}$
is the penalization parameter for the scalar
${\it\xi}$
and
${\it\kappa}$
its diffusivity. In the limit
${\it\eta}_{{\it\xi}}\rightarrow 0$
the solution of the penalized equation tends to the solution of the advection–diffusion equation with no-flux boundary conditions and the penalization error is of
$O(\sqrt{{\it\eta}})$
, as shown in Kadoch et al. (Reference Kadoch, Kolomenskiy, Angot and Schneider2012).
4. Application to fluid turbulence
In the following we consider different applications to two- and three-dimensional hydrodynamic incompressible flows with no-slip walls and also the transport of a passive scalar where no-flux boundary conditions are imposed. The numerical method is based on a pseudo-spectral discretization of the penalized Navier–Stokes and advection–diffusion equations and for details we refer the reader to Schneider (Reference Schneider2005), Kolomenskiy & Schneider (Reference Kolomenskiy and Schneider2009), Kadoch et al. (Reference Kadoch, Kolomenskiy, Angot and Schneider2012) and Engels et al. (Reference Engels, Kolomenskiy, Schneider and Sesterhenn2015).
4.1. Two-dimensional turbulence in a circular container
In Schneider & Farge (Reference Schneider and Farge2005) we presented numerical simulations of two-dimensional decaying turbulence in a circular container with no-slip boundary conditions computed at resolution 10242 using the volume penalization technique. Starting with random initial conditions with Reynolds number
$\mathit{Re}=5\times 10^{4}$
, where
$\mathit{Re}$
is based on the domain size, the turbulent kinetic energy and the kinematic viscosity
${\it\nu}$
, the flow rapidly exhibits self-organization into coherent vortices. Two snapshots of the vorticity field
${\it\omega}$
at later times are shown in figure 3. One-dimensional cuts in figure 4(a) illustrate the intermittent character of the vorticity field and the spikes at the domain boundary show the strong production of vorticity in a thin boundary layer due to the no-slip boundary conditions. The cut of the mask function together with the cuts of the velocity components confirm that in the penalization domain the velocity does indeed vanish. The formation of coherent vortices and the viscous boundary layer have a significant impact on the production and decay of integral quantities. The evolution of kinetic energy
$E(t)$
, enstrophy
$Z(t)$
and palinstrophy
$P(t)$
are shown in figure 4(b). The corresponding balance equations read
$$\begin{eqnarray}d_{t}E=-2{\it\nu}Z,\quad d_{t}Z=-2{\it\nu}P+\oint _{\partial {\it\Omega}}\,{\it\omega}\,(\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\omega})\,\text{d}s,\end{eqnarray}$$
where
$\boldsymbol{n}$
denotes the outer normal with respect to
$\partial {\it\Omega}$
. The source term in the enstrophy dissipation equation involves the vorticity and its gradient at the boundary and yields a significant contribution in the small viscosity limit. The no-slip wall produces vortices which are injected into the bulk flow and tend to compensate the enstrophy dissipation as observed in figure 4. The self-organization of the flow is also reflected by the transition of the initially Gaussian vorticity probability density function (p.d.f.) towards a distribution with exponential tails. Because of the presence of coherent vortices the pressure p.d.f. becomes strongly skewed with exponential tails for negative values. Details can be found in Schneider & Farge (Reference Schneider and Farge2005).
Decaying two-dimensional flow in a container with no-slip boundary conditions at initial Reynolds number
$5\times 10^{4}$
. Snapshots of vorticity
${\it\omega}$
. From Schneider & Farge (Reference Schneider and Farge2005). Reproduced with kind permission of The American Physical Society.
Decaying two-dimensional flow in a container with no-slip boundary conditions at initial Reynolds number
$5\times 10^{4}$
. (a) One-dimensional vertical cuts of vorticity, the two velocity components and the mask function. (b) Evolution of the kinetic energy
$E(t)=(1/2)\int |\boldsymbol{u}|^{2}\,\text{d}\boldsymbol{x}$
, enstrophy
$Z(t)=(1/2)\int |{\it\omega}|^{2}\,\text{d}\boldsymbol{x}$
and palinstrophy
$P(t)=(1/2)\int |\boldsymbol{{\rm\nabla}}{\it\omega}|^{2}\,\text{d}\boldsymbol{x}$
in double logarithmic representation, where
${\it\tau}=t/t_{e}$
is based on the initial eddy turn over time
$t_{e}=\sqrt{2Z(t+0)}=0.061$
. The flow has been integrated for
$650t_{e}$
, corresponding to more than
$10^{5}$
time steps. The inset shows the corresponding log–lin representation. From Schneider & Farge (Reference Schneider and Farge2005). Reproduced with kind permission of The American Physical Society.
4.2. Passive scalar transport in two-dimensional confined turbulence
To illustrate the volume penalization for imposing no-flux boundary conditions we show a numerical simulation of a flow with a passive scalar in a simplified mixing device (Kadoch et al. Reference Kadoch, Kolomenskiy, Angot and Schneider2012). The penalized Navier–Stokes equations (3.12) and (3.13) are solved together with the penalized advection–diffusion equation (3.14) in a square domain with periodic boundary conditions using a Fourier pseudo-spectral method. The fluid domain corresponds to a circular vessel in which a cross-shaped rotor in the centre of the domain rotates in the clockwise direction. The boundary conditions are no-slip for the velocity and no-flux for the passive scalar. The vorticity field in figure 5 illustrates the formation of boundary layers at the wall and at the rotor. The formed vortex sheets destabilize and detach forming coherent vortices which are ejected into the bulk flow. The corresponding passive scalar field, where the initial condition corresponds to a Gaussian blob, is advected by the mean rotation induced by the rotor. The mixing process is further enhanced by the coherent vortices generated by the rotor blades and the boundary layer of the vessel. For further details we refer the reader to Kadoch et al. (Reference Kadoch, Kolomenskiy, Angot and Schneider2012).
(a) Vorticity field in a circular vessel with a cross-shaped clockwise rotating rotor. (b) Corresponding passive scalar field. From Kadoch et al. (Reference Kadoch, Kolomenskiy, Angot and Schneider2012). Reproduced with kind permission of Elsevier Inc.
4.3. Flow past flexible flapping plates in three dimensions
Motivated by simplified models for swimming organisms or robots, which rely on chordwise flexible elastic plates, we present numerical simulations of fluid–structure interaction using the volume penalization. We consider a plate made out of linearly elastic inextensible material, which is perfectly rigid in the spanwise direction but flexible in the chordwise direction. The plate can be thus modelled by the nonlinear beam equation with clamped free boundary conditions and it is solved with classical finite differences. The fluid part is solved with a pseudo-spectral method and volume penalization (Engels et al. Reference Engels, Kolomenskiy, Schneider and Sesterhenn2015). Depending on the fluid/plate density ratio, up to 25 iterations of the fluid–solid coupling are necessary within each time step. Details on the numerical method are described in Engels et al. (Reference Engels, Kolomenskiy, Schneider and Sesterhenn2015).
In Engels et al. (Reference Engels, Kolomenskiy, Schneider and Sesterhenn2014) we considered a configuration with imposed mean flow, and imposed at the leading edge of the plate a sinusoidal pitching motion. The Reynolds number is about
$\mathit{Re}\approx 1000$
. We first simulated a swimmer with a rectangular plate and compared the results with a recent experimental study, before considering also an expanding and a contracting shape of the plate. Flow visualizations for the three geometries showing vorticity isosurfaces are given in figure 6. The tip vortices observed in figure 6 originate from three-dimensional effects due to the finite span. These have important effects for predicting the swimmer’s cruising velocity, since they contribute significantly to the drag force. We found that the cruising velocity of the contracting swimmer is larger than of the rectangular one, which in turn is larger than the expanding one. This observation can be explained by the relative importance of the tip vortices which interact differently with the flexible plates for the three considered geometries of the swimmer. For the contracting case the tip vortices rapidly detach and thus reduce drag, while in the expanding case they are attached down to the trailing edge.
Flow generated by flapping cordwise flexible plates of (a) contracting, (b) rectangular and (c) expanding shape with imposed pitching motion at their leading edge and imposed mean flow (from left to right). Shown are isosurfaces of vorticity
$\Vert {\bf\omega}\Vert =17.5$
. From Engels et al. (Reference Engels, Kolomenskiy, Schneider and Sesterhenn2014).
5. Applications to plasma turbulence
For the numerical simulation of plasma turbulence we consider the non-ideal MHD equations ((3.12)–(3.14)) in which both viscous and resistive effects are taken into account. The magnetic Prandtl number, defined as the ratio between kinematic viscosity
${\it\nu}$
and magnetic diffusivity
${\it\lambda}$
, is equal to one. An isothermal, incompressible plasma is considered with uniform and constant transport coefficients. This approximation simplifies the problem as much as possible, while retaining the required level of complexity to study the nonlinear dynamics. The boundary conditions corresponding to solid domains which are perfect conductors are imposed with the volume penalization method. The numerical code is based on a Fourier pseudo-spectral discretization using FFTs to compute the derivatives and to solve the Poisson equations. It is described in detail in Morales et al. (Reference Morales, Leroy, Bos and Schneider2014b
), including benchmarking and detailed validation studies.
5.1. Spin-up in two-dimensional confined MHD
We consider first two-dimensional decaying MHD turbulence in bounded domains and investigate the spontaneous self-organization with a particular emphasis on the symmetry-breaking induced by the shape of the confining boundaries; for details we refer to Bos et al. (Reference Bos, Neffaa and Schneider2008) and Neffaa et al. (Reference Neffaa, Bos and Schneider2008). This symmetry-breaking is quantified by the angular momentum, which is shown to be generated rapidly and spontaneously from initial conditions free from angular momentum when the geometry lacks axisymmetry. In Bos, Neffaa & Schneider (Reference Bos, Neffaa and Schneider2010) this effect was illustrated by considering circular, square, and elliptical boundaries. It was shown that the generation of angular momentum in non-axisymmetric geometries can be enhanced by increasing the magnetic pressure. Moreover, the effect becomes stronger at higher Reynolds numbers, which are based on the turbulent kinetic energy, the length scale of the domain and the kinematic viscosity. The generation of magnetic angular momentum or angular field, previously observed in Bos et al. (Reference Bos, Neffaa and Schneider2008) at low Reynolds numbers, becomes weaker at larger Reynolds numbers. For axisymmetric geometries, the generation of angular momentum is absent; nevertheless, a weak magnetic field can be observed. The derived evolution equations for both the angular momentum and angular field yield possible explanations for the observed behaviour.
In figure 7 we show simulations of two-dimensional decaying MHD turbulence inside an ovoid and in a D-shaped geometry starting with random initial conditions (Schneider, Neffaa & Bos Reference Schneider, Neffaa and Bos2011). For both geometries we observe that the two-dimensional magnetized plasma self-organizes into a state containing large-scale flow structures, illustrated by the stream function in figure 7 (left) and vorticity, figure 7 (right). To quantify the spin-up we consider the angular momentum defined as
$$\begin{eqnarray}L_{u}(t)=\int _{{\it\Omega}_{f}}\boldsymbol{e}_{\boldsymbol{z}}\boldsymbol{\cdot }\boldsymbol{r}\times \boldsymbol{u}\,\text{d}s=2\int {\it\psi}\,\text{d}s,\end{eqnarray}$$
where
$\boldsymbol{e}_{\boldsymbol{z}}$
is the unit vector in the direction perpendicular to the plane,
$\boldsymbol{r}$
the position vector with respect to the center of the domain and
${\it\psi}$
the stream function.
$L_{u}$
quantifies the fluid rotation. The maximum angular momentum for a given kinetic energy is obtained for a fluid in solid body rotation. In this particular realization we find that the generation of angular momentum is stronger in the ovoid than in the D-shaped geometry, as shown in figure 8.
Two-dimensional magnetized plasma: stream function (left) and vorticity (right) for the ovoid (top) and the D-shaped geometry (bottom). The visualizations correspond to the time instants at which the absolute value of the angular momentum reaches its maximum. From Schneider et al. (Reference Schneider, Neffaa and Bos2011). Reproduced with kind permission of Elsevier Inc.
Two-dimensional magnetized plasma: time evolution of the angular momentum
$L_{u}(t)$
for the ovoid and the D-shaped geometry. From Schneider et al. (Reference Schneider, Neffaa and Bos2011). Reproduced with kind permission of Elsevier Inc.
Confined MHD flows in toroidal domains. Asymmetric poloidal geometry. Streamlines coloured by toroidal velocity value for
$M=7.5$
(a) and
$M=75.2$
(b). From Morales et al. (Reference Morales, Bos, Schneider and Montgomery2012). Reproduced with kind permission of The American Physical Society.
(a) Sketch of the toroidal geometry with asymmetric poloidal cross-section, (b) toroidal velocity field component isocontours (blue
$+9\times 10^{-4}$
, orange − 9 × 10−4) and (c) perturbed toroidal magnetic field component isocontours (red
$+0.025$
, orange
$+0.04$
, yellow
$+0.05$
). From Morales et al. (Reference Morales, Leroy, Bos and Schneider2014b
). Reproduced with kind permission of Elsevier Inc.
Confined MHD flows in toroidal domains. (a) Normalized toroidal velocity
$\langle u_{{\it\theta}}^{2}\rangle /\langle \boldsymbol{u}^{2}\rangle$
. The inset shows the modulus of the cosine of the angle between the velocity and the magnetic field. (b) Volume-averaged toroidal angular momentum
$\langle L_{{\it\theta}}\rangle$
. All quantities are plotted as a function of the viscous Lundquist number
$M$
for the symmetric and asymmetric poloidal cross-section geometries. From Morales et al. (Reference Morales, Bos, Schneider and Montgomery2012). Reproduced with kind permission of The American Physical Society.
5.2. Self-organization of confined MHD flows in toroidal domains
The spatio-temporal self-organization of visco-resistive magnetohydrodynamics in a toroidal geometry (see figure 10
a below) was studied in Morales et al. (Reference Morales, Bos, Schneider and Montgomery2012) using fully three-dimensional simulations considering two geometries: a torus with a symmetric poloidal cross-section; and one with an asymmetric poloidal cross-section. The magnetized plasma is initially in a quiescent state and curl-free toroidal magnetic and electric fields are imposed. The simulations show that spontaneously a flow field is generated in both geometries and the magnetized plasma starts to move. Moreover, the flow evolves from dominantly poloidal to toroidal when the Lundquist (or Reynolds) numbers
$M$
are increased. Here the viscous Lundquist number
$M$
is defined as
$M=C_{A}L/{\it\nu}$
, where
$C_{A}$
is the toroidal Alfvén velocity,
$L$
the diameter of the cross-section of the torus and
${\it\nu}$
the kinematic viscosity. In Morales et al. (Reference Morales, Bos, Schneider and Montgomery2012) we have shown that this toroidal organization of the flow is consistent with the tendency of the velocity field to align with the magnetic field. Furthermore, we found that the up–down asymmetry of the geometry causes the generation of a non-zero toroidal angular momentum.
Figure 9 illustrates the streamlines, coloured by the toroidal velocity value for two Lundquist numbers for the asymmetric poloidal geometry. For the low Lundquist number case (figure 9 a) we do indeed observe a pair of counter-rotating vortices in the poloidal plane, while for the larger Lundquist number (figure 9 b) the flow starts moving in the toroidal direction. A similar behaviour is observed for the symmetric poloidal cross-section.
The velocity and the perturbed magnetic toroidal component isocontours at the steady state for a toroidal geometry with asymmetric poloidal cross-section are shown in figure 10(b,c). The perturbed toroidal magnetic field is created by the velocities in the poloidal plane. This component of the magnetic field is important because it generates a toroidal Lorentz force that induces the toroidal velocities.
Figure 11(a) shows that the toroidal velocity increases with the viscous Lundquist number in both geometries and saturates at about 86 % of the total squared speed. The inset shows the modulus of the cosine of the angle between the velocity and the magnetic field and thus quantifies that the velocity fields tend to align with the magnetic field for increasing
$M$
. However, for the volume-averaged toroidal angular momentum, defined as
$\langle L_{{\it\theta}}\rangle =(1/V)\int _{V}Ru_{{\it\theta}}\,\text{d}V$
, we observe fundamental differences for the two geometries as shown in figure 11(b). For the torus with the symmetric poloidal cross-section,
$\langle L_{{\it\theta}}\rangle$
identically vanishes for all considered Lundquist numbers. In contrast we observe for the asymmetric case that the toiroidal angular momentum increases with the Lundquist number.
5.3. Effect of toroidicity in RFP dynamics
Finally, we consider the reversed field pinch (RFP) dynamics and study the influence of the curvature of the imposed magnetic field. In RFP experiments the plasma evolved to quasi-stationary equilibria characterized by Beltrami minimum energy states for which the magnetic field corresponds to eigenfunctions of the curl operator (Nordlund Reference Nordlund1994). In Morales et al. (Reference Morales, Bos, Schneider and Montgomery2014a
) we compared the RFP flow of a magnetofluid in a torus with aspect ratio 1.83 with the flow in a periodic cylinder, and studied the persistence of these dominant helical modes for varying pinch ratios. The flow configuration of both geometries is illustrated in figure 12. The pinch ratio is defined as the wall-averaged poloidal magnetic field divided by the volume-averaged toroidal magnetic field,
${\it\Theta}=\overline{B_{P}}/\langle B_{T}\rangle$
. The ratios of kinetic energy of the dominant axial and toroidal mode of the total kinetic energy versus the pinch ratio are shown in figures 13(a) and 13(b) for the cylindrical and the toroidal geometry, respectively. We find that an axisymmetric toroidal mode is always present in the toroidal, but absent in the cylindrical configuration. In particular, in contrast to the cylinder, the toroidal case presents a double poloidal recirculation cell with a shear localized at the plasma edge. Quasi-single-helicity states are found to be more persistent in toroidal than in periodic cylinder geometry. The dominant helical modes at pinch ratios
${\it\Theta}\geqslant 2$
are illustrated in figure 13 by showing axial and toroidal velocity isosurfaces to provide further insight into the flow topology. For further details we refer to Morales et al. (Reference Morales, Bos, Schneider and Montgomery2014a
).
Flow configuration of the RFP dynamics: (a) toroidal geometry; (b) periodic cylinder. From Morales et al. (Reference Morales, Bos, Schneider and Montgomery2014a ). Reproduced with kind permission of IOP Publishing Ltd.
RFP dynamics. Ratio of the kinetic energy of (a) the dominant axial and (b) toroidal modes over the total kinetic energy for (a) the cylindrical and (b) the torus geometry for
$M=888$
as a function of the pinch ratio
${\it\theta}$
. Also shown are visualizations of the modes: (a) axial velocity isosurfaces
$+0.008$
(blue) and
$-0.008$
(orange); (b) toroidal velocity isosurfaces
$+0.007$
(blue) and
$-0.007$
(orange). From Morales et al. (Reference Morales, Bos, Schneider and Montgomery2014a
). Reproduced with kind permission of IOP Publishing Ltd.
6. Conclusion and perspectives
We have reviewed immersed boundary methods with a special focus on the volume penalization method for imposing Dirichlet (corresponding to no-penetration and no-slip conditions) or Neumann boundary conditions (corresponding to no-flux conditions) in complex geometries. The mathematical concepts for choosing the parameters involved, i.e., the penalization parameter and the grid size, have been illustrated by considering simple one-dimensional problems. Applications to hydrodynamic and magnetohydrodynamic turbulence in complex geometries using a simple Fourier pseudo-spectral method, which can be parallelized on massively parallel machines using standard libraries like P3DFFT, illustrated the versatile use of this technique for various problems encountered in computational physics. Toroidal geometries can thus be efficiently handled, even including asymmetric poloidal cross-sections, and simulations for higher Lundquist numbers become feasible. An essential feature of the volume penalization method is that it becomes more attractive for computing fluid flows for small viscosity values, i.e., for high Reynolds/Lundquist number flows. The reason is that the effective penalization boundary layer size depends on the product of viscosity and the penalization parameter and thus the method becomes more precise without using prohibitively small penalization parameters, cf. Nguyen van yen et al. (Reference Nguyen van yen, Kolomenskiy and Schneider2014).
One perspective of current research is the development of higher-order penalization methods which allow faster convergence to be obtained. The Cartesian grid introduces in dimension larger equal to two a staircase effect for complex (non-grid aligned) geometries and the approximation of the mask function thus reduces to first order. Techniques based on interpolation to obtain higher order for complex geometries have been proposed, e.g., in Sarthou et al. (Reference Sarthou, Vincent, Caltagirone and Angot2008) in the context of finite volume formulations. Another challenging topic are more sophisticated boundary conditions for MHD flows, overcoming the limitation of perfect conductors in the surrounding solid domain.
Acknowledgements
K.S. is grateful to S. Benkadda for inviting him to give a review lecture at the ITER International School 2014 ‘High performance computing in fusion science’ on which the manuscript is based on. K.S. is also indebted to P. Angot, W. Bos, T. Engels, M. Farge, B. Kadoch, D. Kolomenskiy, M. Leroy, J. Morales, S. Neffaa, R. Nguyen van yen, M. Roberts with whom the volume penalization method has been developed together and published in numerous research papers on which this review is based on. K.S. gratefully acknowledges financial support from the ANR project SiCoMHD (ANR-Blanc 2011-045) and thanks Matteo Faganello for useful comments on the manuscript.
