Comparison between direct numerical simulations and effective models for fluid-porous flows using penalization

Abstract

This work is devoted to a numerical study of two-dimensional incompressible flows in a fluid-porous medium system placed on an impermeable wall. Flow in the whole system is studied either at the pore scale or at the Darcy scale (using a one-domain approach) and the models at both scales are solved with a penalization method using the same formal Navier–Stokes equations modified by a Darcy-like term. Several effective medium penalized models are considered for the simulations. Various flow regimes are investigated ranging from laminar to turbulent. The velocity profiles inside and outside the porous medium show significant discrepancies between the different penalization models compared to the direct numerical simulations. This work motivates further studies about the spatial variations of the penalization coefficient to be introduced in the effective medium models in order to better reproduce the physics near the fluid-porous medium boundaries.

Appendices

Appendix 1

The values of the components of the tensor ${{\mathbf {\mathsf{{H}}}}}^{-1}$ are computed for a square pattern of parallel cylinders of square cross section having a porosity $\epsilon =0.913$ and a flow orthogonal to the cylinders axes. The components $H_{xx}$, $H_{yx}$, $H_{xy}$ and $H_{yy}$, of tensor ${{\mathbf {\mathsf{{H}}}}}$ are satisfying the following properties:

$$\begin{aligned} H_{xx}(\theta ,Re_{local}) &= H_{yy}(90^\circ -\theta ,Re_{local}),\quad \; H_{yx}(\theta ,Re_{local}) = H_{xy}(90^\circ -\theta ,Re_{local}) \\ H_{xy}(\theta ,Re_{local})&= H_{yx}(90^\circ -\theta ,Re_{local}),\quad \; H_{yy}(\theta ,Re_{local}) = H_{xx}(90^\circ -\theta ,Re_{local}). \end{aligned}$$

Then, for $\theta$ from $95^\circ$ to $180^\circ$ the components of the tensor are given by:

$$\begin{aligned} H_{xx}(\theta ,Re_{local})&= H_{xx}(180^\circ -\theta ,Re_{local}),\quad \; H_{yx}(\theta ,Re_{local}) = -H_{yx}(180^\circ -\theta ,Re_{local}) \\ H_{xy}(\theta ,Re_{local})&= -H_{xy}(180^\circ -\theta ,Re_{local}),\quad \; H_{yy}(\theta ,Re_{local}) = H_{yy}(180^\circ -\theta ,Re_{local}). \end{aligned}$$

Finally, for $\theta$ from $185^\circ$ to $360^\circ$ the components of the tensor are given by:

$$\begin{aligned} H_{xx}(\theta ,Re_{local})&= H_{xx}(360^\circ -\theta ,Re_{local}),\quad \; H_{yx}(\theta ,Re_{local}) = -H_{yx}(360^\circ -\theta ,Re_{local}) \\ H_{xy}(\theta ,Re_{local})&= -H_{xy}(360^\circ -\theta ,Re_{local}), \; H_{yy}(\theta ,Re_{local}) = H_{yy}(360^\circ -\theta ,Re_{local}). \end{aligned}$$

The same symmetries are applicable to the tensor ${{\mathbf {\mathsf{{H}}}}}^{-1}$, the component values of which are provided as additional material.

Appendix 2

In this appendix, a heuristic modification of the K model is explored. The idea is to investigate the possibility of keeping the same form of the effective-medium momentum equation, which has the advantage of being rather simple and solved with the same numerical code, by empirically modifying the Darcy-like penalization term.

An examination of the flow around the front of the porous region with the K model reveals that there is almost no transition between the fluid and the porous medium (see Fig. 10 left) at $Re=1000$. The zoom of the velocity field corresponds to the same part of the domain shown in Fig. 3. Indeed, for the K model the mean velocity inside the porous domain (averaged with respect to y, which varies from 0 to 30) near the porous medium entrance is quite large. It ranges from 0.83 at $x=100$ to 0.14 at $x=200$ (Table 3) instead of 0.33 at $x=100$ and 0.01 at $x=200$ for the reference flow. Data reported in Table 3 show that due to the re-circulation zone on top of the porous medium, the mean velocity inside the porous region becomes negative in the second part of this re-circulation zone, namely between $x=250$ and $x=300$. However, the velocity is always positive with the K model.

Fig. 10

Zoom at the same location of the instantaneous velocity field around the front of the porous rectangle for $Re=1000$. Left: the K model, right: the improved K–U model. Both models are computed on the G6 grid. The length and color of the velocity vectors are related to the velocity magnitude from 0 (dark blue) to the maximum (dark red). (Color figure online)

In this model, the coefficient $k/\epsilon$ in ${\varvec{\Gamma }}$ is equal to 0.0453 in the present computation. For $Re=1000$, it yields ${\varvec{\Gamma }}^{-1}=45.3{{\mathbf {\mathsf{{I}}}}}$. This value seems correct in the middle part of the porous medium but does not perform well near the fluid-porous medium boundary, specially when there is a strong normal velocity ahead. To decrease the velocity inside the porous medium close to the surface, it is necessary to decrease locally the value of ${\varvec{\Gamma }}^{-1}$. This can be done, for example, by arbitrarily dividing ${\varvec{\Gamma }}^{-1}$ by an order of magnitude at the boundary and then increasing it linearly in the first $10\%$ of the rectangle size to reach the right value ${\varvec{\Gamma }}^{-1}=45.3{{\mathbf {\mathsf{{I}}}}}$ and keeping it until the exit section of the porous medium. This can be done by replacing ${\varvec{\Gamma }}^{-1}$ by ${\tilde{\varvec{\Gamma }}}^{-1}$ defined as ${\tilde{\varvec{\Gamma }}}^{-1}=\frac{4}{144-x}{\varvec{\Gamma }}^{-1}$ for $100\le x\le 140$. The resulting model is called the K–K model in the following.

Table 3 Comparison of the mean velocity inside the porous rectangle between the reference flow computed by DNS on the $\widetilde{G8}$ grid and the three K models computed on the G6 grid at various locations

The solution computed with the K–K model has a mean value very close to the reference flow in the entrance region of the porous medium, mostly for $100\le x\le 200$. This is very encouraging but the flow computed with the K–K model is still quite different from the reference flow further into the porous medium (see Fig. 11). An alternative way to modify the penalization coefficient for the effective medium K model is to relate this coefficient to the local velocity magnitude. Then the penalization coefficient can be modulated not only in the front part of the porous medium but within the whole porous medium. With this approach it is possible, in principle, to also encompass the situation when a strong vortex travels down on top of the porous medium. When the velocity magnitude exceeds a chosen threshold, the penalization coefficient is reduced accordingly. Analyzing the velocity inside the porous medium, this can be applied when the velocity magnitude exceeds a value around 0.1. Therefore, a modified value of ${\varvec{\Gamma }}$ can be proposed as ${\tilde{\varvec{\Gamma }}}^{-1}=\frac{1}{100(||{\mathbf {u}}||-0.1)}{\varvec{\Gamma }}^{-1}$ if $||{\mathbf {u}}||>0.11$. The resulting model is called the K–U model in the following.

Fig. 11

Mean velocity profiles at three different vertical sections positioned along the x-direction at 200, 400 and 500 for $Re=1000$. Left: zoom around the front of the porous rectangle, right: zoom in the middle part above the porous rectangle. The two penalization models computed on the G6 grid are compared to the reference flow computed by DNS on the $\widetilde{G8}$ grid

As shown in Table 3, the K–U model does not perform better than the K–K model in the front part of the porous medium. However, the predictions are much improved downstream as the K–K model yields a very poor approximation (Fig. 11 left). Moreover, the approximation is much better with the K–U model in the whole domain at the exit section (Fig. 11 right). The empirical improvements made on the K model, make possible to find a coherent behaviour at the front of the porous body as shown in Fig. 10 right. The porous rectangle is now clearly visible as well as the re-circulation above like in Fig. 3 right for the DNS. As expected the two snapshots can not be directly superimpose because they represent instantaneous solutions that are necessarily different with two approaches so far apart which have a different transient. In addition the Kelvin-Helmholtz vortices are smaller on the finest $\widetilde{G8}$ grid. However, what is captured by the effective models is the mean flow behaviour.

Fig. 12

Mean velocity profiles at three different vertical sections positioned along the x-direction at 200, 400 and 500 for $Re=1000$. Left: zoom around the front of the porous rectangle, right: zoom in the middle part above the porous rectangle. The K–U model computed on the G6 grid is compared to the reference flow computed by DNS on the $\widetilde{G8}$ grid and to the H model computed on the G7 grid

To conclude this section, it is worth comparing the performance of the K–U model computed on the G6 grid and the best one obtained so far, i.e., the H model computed on the G7 grid, with respect to the DNS reference obtained on the $\widetilde{G8}$ grid for $Re=1000$. In Fig. 12, the corresponding horizontal velocity profiles, u, are presented. In the porous region (left graph of Fig. 12), the modification on ${\varvec{\Gamma }}$ in the K–U model improves significantly the prediction of u at the entrance, although a re-circulation is obtained which does not correspond to the reference flow. The improvement is less significant at the center of the porous region whereas the prediction deteriorates at the exit. Above the porous region, (right part of Fig. 12), the K–U model improves the prediction of u in all the three investigated sections, in particular at $x=200$. The effective integration of the interface, in the case of high speed, to compute the full tensor ${{\mathbf {\mathsf{{H}}}}}$ should provide results close to the DNS. Finally, it must be pointed out that the K–U model only requires $2\%$ of the DNS CPU time.