NUMERICAL SIMULATION OF WALL-BOUNDED FLOWS USING A SPECTRAL METHOD WITH VOLUME PENALIZATION

The volume penalization method, which allows to impose no-slip boundary conditions, is assessed for wall-bounded flows. For the numerical solution of the penalized equations a spectral method is used. Considering a two-dimensional Poiseuille flow, the solution of the Navier-Stokes penalized equation is computed analytically and the convergence of the numerical solution is studied. To illustrate the properties of the approach we compute a three-dimensional turbulent channel flow imposing a constant flow rate. The obtained results are compared with reference data of Kim et al. [10]. Résumé. La méthode de pénalisation en volume, qui permet d’imposer des conditions aux limites non glissement, est analysée pour des écoulements avec parois. Pour la solution numérique des équations pénalisées une méthode spectrale est utilisée. Nous considerons un écoulement bidimensionnel de type Poiseuille et la solution exacte des équations pénaliées est calculée analytiquement. Ensuite nous étudierons la convergence de la solution numérique. Des simulations d’un écoulement tridimensionnel turbulent dans un canal avec un débit constant montrent les propriétés de cette approche.


Introduction
Efficient and accurate numerical simulation of turbulent flows in complex and even time varying geometries, is still a major challenge in computational fluid dynamics.Immersed boundary methods, and volume penalization methods in particular, have obtained growing interest during the last decades in this context.For reviews, we refer to Mittal & Iaccarino [1], Peskin [2] and Schneider [3].Schneider [3] also illustrates some applications for magneto-hydrodynamic and plasma turbulence.Instead of solving the governing equations, typically Navier-Stokes equations of either compressible or incompressible fluid, in the complex geometry using body fitted grids or coordinate transforms, the idea is to solve the equations in a simple, e.g.Cartesian geometry.The boundary conditions are then imposed by appending additional terms to the equations.In most of the preceding applications flows, primarily flows in the low Reynolds number regime have been considered.Some exceptions can be found for two-dimensional (2D) confined turbulence in Schneider & Farge [4] and recently also for three-dimensional (3D) turbulent flow over periodic hills in Kadoch et al. [5].
The volume penalization method is physically motivated, as the wall, or solid obstacles are modeled by a porous medium.Inside the porous medium the flow is governed by Darcy's law.In the limit when the permeability of the porous medium tends to zero, we recover a solid wall or a solid obstacle and the enforced boundary conditions correspond indeed to no-slip conditions.This has been shown rigorously for Navier-Stokes equations in Angot et al. [6] and Carbou & Fabrie [7].Further analyses of the penalized problems can be found in Nguyen van yen et al. [8].Detailed numerical validation tests of the volume penalization method in the context pseudo-spectral discretization and for simple and complex geometries can be found in Refs.[8,9].
The verification of the volume penalization for high Reynolds number is however still in its infant state and the current work can be seen as a contribution into this direction.Here we consider moderate Reynolds number flows under a mean pressure gradient and we perform analysis and numerical simulation for two canonical flows.For a 2D Poiseuille flow we derive the exact solution of the penalized equations and perform estimations of the numerical errors.For a 3D turbulent channel flow we propose a method for keeping the flow rate constant, and also a technique for determining the frictional velocity at the walls.The turbulence statistics are also assessed and compared with available DNS data from Kim et al. [10] computed with a Fourier-Chebyshev spectral method.
The paper is organized as follows: Section 2 presents the basic equations using the penalization method.In Sections 3 and 4, we consider two canonical flows: one is the plain 2D Poiseuille flow and the other is the 3D turbulent channel flow.A detailed numerical verification is presented.Conclusions and perspectives are given in Sec. 5.

Basic equations using penalization method
We consider motions of the incompressible fluid of unit density subjected to a mean pressure gradient and bounded by solid walls in the Cartesian coordinate system x = (x, y, z).The walls are modeled by the use of the volume penalization method.The solid and fluid regions are respectively denoted by Ω s and Ω f (see Fig. 2.1).The governing equations of the motions are given as follows: where u(x, t) = (u, v, w) is the velocity, p(x, t) is the pressure fluctuation, G is the intensity of the mean pressure gradient in the x-direction, e x = (1, 0, 0), ν is the kinematic viscosity, t is time and ∇ = (∂/∂x, ∂/∂y, ∂/∂z).
The penalization parameter η denotes the permeability of the solid region Ω s .The mask function χ(x, t), describing the geometry of Ω s , is defined as where Ω s denotes the closure of Ω s , and the whole region Ω is defined by Ω = Ω f ∪ Ω s .Here and hereafter, we omit the arguments x and t, unless otherwise stated.

Plain Poiseuille flow using penalization
In this section, we consider 2D plain Poiseuille flows using the volume penalization method.The two walls are parallel to the x-direction and normal to the y-direction.Thus, the solid domain Ω s is defined by Ω s = {(x, y)|0 ≤ x < 4π, 1 ≤ |y| < L}, while the fluid domain Ω f is defined by Ω f = {(x, y)|0 ≤ x < 4π, −1 < y < 1}.We employ 2L-periodic boundary condition in the y-direction and 4π-periodic boundary condition in the xdirection.
We first obtain an exact solution of the steady parallel flow subjected to the mean pressure gradient, which satisfies the 2D representation of Eqs. ( 1) and (2).Then we perform numerical simulations with L = 1.28 and assess the numerical error in the penalization method.

The exact solution of the penalized equations
In the following, we compute the exact solution of the penalized equations and show the convergence towards the solution of the non-penalized equation in the limit when the penalization parameter tends to zero.The steady parallel flow is given by u = (u(y), 0) and p = 0.It thus satisfies Imposing C 1 -continuity conditions at y = ±1, and using the 2L-periodic boundary condition in the y-direction, we obtain after some algebra the exact solution of the penalized equations, where . The parameter λ denotes the thickness of the penalization boundary layer.For η → 0, we find which shows that in the limit of vanishing permeability, i.e. for the solid walls, we obtain indeed the 2D Poiseuille flow.It can be noted that u(±L) → 0 as L → ∞.

Numerical verification
Now we simulate 2D steady flow using the volume penalization method for L = 1.28, ν = 0.1 and G = 2.0.To study the numerical convergence, we use different number of grid points in the y-direction, N y , and different penalization parameters η.The number of grid points in the x-direction is fixed, i.e.N x = 32.The corresponding Reynolds number defined by G/(2ν 2 ) is 100.We use a Fourier pseudo-spectral method where the aliasing errors are removed by the 2/3-rule.The first-order explicit Euler method is employed for time marching.We regard the flows as steady ones, when the maximum value of |∂ û(k)/∂t| + |∂v(k)/∂t| in terms of the wave vector k becomes less than 10 −6 .Here, u(k) and v(k) are the Fourier transforms of u and v, respectively.The errors between the numerical solutions and the exact solution, Eq. ( 7), are measured by the L 2 and L ∞ norms.Figures 3.1 and 3.2 present the N y -and η-dependence of the errors, respectively.Figure 3.1 shows first-order convergence of the errors in terms of N y .In Fig. 3.2, we observe that the errors decay approximately as η 1/2 for η 10 −3 , while for smaller η, they are found to be independent of η.

Numerical settings and penalized representation of physical quantities
We move on to the numerical simulation of turbulent channel flows using the volume penalization method.The computational domain is periodic in each Cartesian direction and it consists of the fluid region Ω f and the (7), in the whole region at Ny = 256.The green and red lines denote the 2 and ∞ norms, respectively.The blue dotted line denotes η 1/2 decay.solid region Ω s , which are respectively defined as where L x = 4π, L y = 2.56 and L z = 4π/3.The x-, y-and z-directions are respectively streamwise, wall-normal and spanwise directions.Again we use a fully dealiased Fourier pseudo-spectral method, now with a second-order Runge-Kutta method for time marching.The aliasing errors are removed by the 2/3 rule.In order to compare our numerical results with a reference direct numerical simulation (DNS) results by Kim et al. [10], we chose the kinematic viscosity ν such that the Reynolds number based on the mean velocity over the fluid domain, defined by Re m = 2U m /ν, yields 5600.The mean velocity is defined by U m = 1 −1 u (y)dy/2, where • denotes the y-dependent spatial average of • over the x-z plane.The total number of the grid points is N x N y N z , where N α (for α = x, y, z) is the number of grid points in the α-direction.We perform all simulations with N x = N z = 128, but to study the convergence we use different N y and different values of the penalization parameter η.A summary of the parameters can be found in Table 1.The number of grid points in the y-direction is set to N y = 512, 1024, or 2048 and the corresponding grid width is L y /N y .Note that the grid with N y = 1024 is similar to the minimum grid width in the Chebyshev reference simulations, i.e., 8π 2 /N 2 KMM , where N KMM = 129 is the number of grid points in the Chebyshev case in Ref. [10].The value of the mean pressure gradient term G is determined such that the total flow rate is kept constant.Using the Fourier transformation of Eqs. ( 1) and (2), and taking into account that u(0) is the total flow rate, we obtain We analyze 20 snapshots of data for the turbulent channel flow with intervals of 0.1 washout times.Then we take the ensemble-average over those 20 snapshots to guarantee well-converged statistical results.The averaging starts after the energy has become quasi-stationary.The washout time is defined by 4π/U m .The friction Reynolds number Re τ , defined by Re τ = u τ /ν, is about 180, which is almost the same as Kim et al. [10]   The frictional velocity is given by where u are the streamwise velocity fluctuations defined by u = u − u .The friction velocity u τ is obtained under the assumption that the flows are statistically stationary.After taking the y-dependent spatial average over x − z planes, the x-component of Eq. ( 1) results in Here, we have used uv = u v because v = 0.Then, integrating Eq. ( 13) with respect to y over the range −L y /2 ≤ y ≤ 0, yields Eq. ( 12).We have imposed the conditions that | u v (−L y /2)| 1 and |∂ u /∂ y | 1.These conditions show that all terms related to the velocity fluctuations are sufficiently small at y = −L y /2, if we take an appropriate solid region depending on L y and η.Quantities with the superscript + are expressed in the wall units, i.e., they are non-dimensionalized by u τ and ν.

Numerical results
Visualizations of the isosurface values of the modulus of vorticity at a given time instant, for the flows computed with different wall normal resolutions, are shown in Fig. 4.1.We observe that both flows exhibit intense vortex tubes near the walls, as observed in previous DNS (e.g., Ref. [11]).Here, y + = u τ (y + 1)/ν.As N y increases, we find that u shows better agreement with the reference data by Kim et al. [10] in the viscous layer, as well as in the core region.The mean pressure p also better agrees with the reference with increasing N y , though for the lowest resolution N y = 512, p slightly oscillates.
In Fig. 4.3, we see the mean distributions, u and p , vs. y + for different η and for N y = 2048.As η becomes smaller, both distributions exhibit better agreement with the reference data.For the largest η = 10 −1 , we observe some discrepancy with the reference data, in particular near the wall, i.e. for y + 10.Now, we examine the root-mean-square (RMS) of the velocity and pressure fluctuations, denoted by u + RMS , v + RMS , w + RMS and p + RMS .Their fluctuations are given by u− u and p− p .Figures 4.4 and 4.5 show respectively the N y -dependence and the η-dependence of these RMS values as a function of y.We observe that all velocity RMSs excellently agree with the reference DNS data [10] near the walls.In contrast, in the core region, these velocity RMS values, u + RMS , v + RMS and w + RMS , are substantially smaller than the reference RMS ones, irrespective of η and N y studied here, though the former velocity RMS exhibit similar y-dependence as the reference ones.These smaller RMS could be attributed to the finite values of the permeability η.In Ref. [12] it was shown that the existence of a porous wall reduces indeed turbulent energy production.Concerning the pressure RMS p + RMS , we find as well some departure from the reference data for all y.

Conclusions
The volume penalization method in combination with a spectral method has been assessed for wall-bounded flows in simple geometries.The no-slip boundary conditions are modeled by extending the computational domain.Inside the fluid domain the classical incompressible Navier-Stokes equations are solved, while in the added porous domain, which models the walls, the Darcy equation is solved.In the limit when the permeability parameter of the porous domain tends to zero, we recover the Navier-Stokes solution with no-slip boundary conditions.
To perform some analysis we considered a 2D Poiseuille flow first.An exact solution of the penalized equations has been computed and the penalization error was shown to be of order η 1/2 .Then we studied the convergence of the numerical solution and found first order convergence in space towards the 2D Poiseuille flow.
Secondly we applied the method to a 3D turbulent channel flow where we imposed a constant flow rate.We analyzed the turbulent statistics and compared the results with the reference DNS data of [10].We found that mean velocity and pressure profiles show better agreement with the DNS results, either with increasing grid points in the wall normal directions, or with decreasing the permeability of the solid region.The RMS values of velocity fluctuations showed similar distributions to the reference data.However, we observed some reduction of their values.It can be speculated that the finite permeability is responsible for reducing the turbulent energy production.The RMS values of the pressure fluctuations exhibit also some departure from the reference computation.
In conclusion the volume penalization is an attractive method which can be easily implemented into existing codes using, e.g.Cartesian grids and classical discretizations.The flow geometry is completely described by a mask function and even for complex and time varying geometries no complicated gridding or mesh generation becomes necessary.Thus the method is very flexible and useful for many applications, e.g. for biological flows encountered in insect flight as shown in [9] or for magnetohydrodynamics in toroidal geometries as presented in [3].Detailed validation tests can be found in the respective papers [3,8,9] where likewise pseudo-spectral discretizations have been used.
A drawback of the penalization method is the reduced convergence order of the method (formally first order only), which is due to the numerical boundary layer.Adaptive discretizations can alleviate this problem, as the grid can be automatically refined near the boundary, as shown by Okamoto et al. [13].The time step restriction, which is physically imposed and not only by the stability requirements of explicit discretizations, cannot be overcome by implicit time schemes.However, for turbulent flows, for which the viscosity is small, the penalization parameter is not required to assume very small values and hence the time step is typically imposed by the CFL condition.
Perspectives of the current work are the extension to multiphysics problems, like coupled problems with thermal diffusion in complex geometries, which may even be time dependent.

Figure 3 . 1 .
Figure 3.1.Ny-dependence of the error between the numerical solution and the exact solution, Eq. (7), in the whole region Ω = Ω f ∪ Ωs at η = 10 −2 .The green and red lines denote the 2 and ∞ norms, respectively.The blue dotted line denotes the N −1 y decay.

Figure 3
Figure 3.2.η-dependence of the error between the numerical solution and the exact solution, Eq. (7), in the whole region at Ny = 256.The green and red lines denote the 2 and ∞ norms, respectively.The blue dotted line denotes η 1/2 decay.

Figure 4 . 2 .Figure 4 . 3 .
Figure 4.2.The Ny-dependence of the mean streamwise velocity u and the mean pressure p at η = 10 −2 .Red, green and blue lines denote the results for Ny = 512, 1024 and 2048, respectively, while the purple points correspond to the reference data by Kim et al. [10].

Figure 4 .
Figure 4.2 shows the y + -dependences of the mean streamwise velocity u and pressure p for different N y with η = 10 −2 .Here, y + = u τ (y + 1)/ν.As N y increases, we find that u shows better agreement with the reference data by Kim et al.[10] in the viscous layer, as well as in the core region.The mean pressure p also better agrees with the reference with increasing N y , though for the lowest resolution N y = 512, p slightly oscillates.In Fig.4.3, we see the mean distributions, u and p , vs. y + for different η and for N y = 2048.As η becomes smaller, both distributions exhibit better agreement with the reference data.For the largest η = 10 −1 , we observe some discrepancy with the reference data, in particular near the wall, i.e. for y + 10.Now, we examine the root-mean-square (RMS) of the velocity and pressure fluctuations, denoted by u + RMS , v + RMS , w + RMS and p + RMS .Their fluctuations are given by u− u and p− p .Figures 4.4 and 4.5 show respectively the N y -dependence and the η-dependence of these RMS values as a function of y.We observe that all velocity RMSs excellently agree with the reference DNS data[10] near the walls.In contrast, in the core region, these velocity RMS values, u + RMS , v + RMS and w + RMS , are substantially smaller than the reference RMS ones, irrespective of η and N y studied here, though the former velocity RMS exhibit similar y-dependence as the reference ones.These smaller RMS could be attributed to the finite values of the permeability η.In Ref.[12] it was shown

Figure 4 . 4 .
Figure 4.4.The Ny-dependence of RMS values of velocity components, u + RMS , v + RMS , w + RMS , and the values of pressure p + RMS for η = 10 −2 together with the reference DNS results of Ref. [10].Red, green and blue lines denote the results for Ny = 512, 1024 and 2048, respectively.The purple points correspond to the reference data.

Table 1 .
Numerical parameters for the different computations of the turbulent channel flow.