Homogenized Flux-Body Force Treatment of Compressible Viscous Porous Wall Boundary Conditions

TL;DR

This paper introduces a homogenization method for modeling porous wall boundary conditions in compressible viscous flow simulations, eliminating the need for pore-resolved meshes and heuristic parameters, and validated through supersonic flow tests.

Contribution

It presents a novel homogenization approach that models porous walls without pore-resolved meshes or heuristic discharge coefficients, including a body force term for friction effects.

Findings

Method accurately predicts supersonic flow over porous walls.

Validated against pore-resolved numerical simulations.

Applicable to practical flow simulations with porous boundaries.

Abstract

A homogenization approach is proposed for the treatment of porous wall boundary conditions in the computation of compressible viscous flows. Like any other homogenization approach, it eliminates the need for pore-resolved fluid meshes and therefore enables practical flow simulations in computational fluid domains with porous wall boundaries. Unlike alternative approaches however, it does not require prescribing a mass flow rate and does not introduce in the computational model a heuristic discharge coefficient. Instead, it models the inviscid flux through a porous wall surrounded by the flow as a weighted average of the inviscid flux at an impermeable surface and that through pores. It also introduces a body force term in the governing equations to account for friction loss along the pore boundaries. The source term depends on the thickness of the porous wall and the concept of an equivalent single pore. The feasibility of the latter concept is demonstrated using low-speed permeability test data for the fabric of the Mars Science Laboratory parachute canopy. The overall homogenization approach is illustrated with a series of supersonic flow computations through the same fabric and verified using supersonic, pore-resolved numerical simulations.

Tables

Table 1. Table 1 : Supersonic porous flow model problem: boundary and initial conditions (additionally, v = w = 0 ​ m ​ s − 1 𝑣 𝑤 0 m superscript s 1 v=w=0\ \mathrm{m\ s^{-1}} in all cases and both subdomains).

$V_s$	Upstream boundary/initial conditions	Downstream boundary/initial conditions
$279.4m{s}^{-1}$	$\rho =0.0076\mathrm{kg}{m}^{-3}$, $u=127.98m{s}^{-1}$, $p=260\mathrm{Pa}$	$\rho =0.004119\mathrm{kg}{m}^{-3}$, $u=0m{s}^{-1}$, $p=112.73\mathrm{Pa}$
$335.9m{s}^{-1}$	$\rho =0.0076\mathrm{kg}{m}^{-3}$, $u=213.31m{s}^{-1}$, $p=260\mathrm{Pa}$	$\rho =0.002774\mathrm{kg}{m}^{-3}$, $u=0m{s}^{-1}$, $p=61.240\mathrm{Pa}$
$367.1m{s}^{-1}$	$\rho =0.0076\mathrm{kg}{m}^{-3}$, $u=255.97m{s}^{-1}$, $p=260\mathrm{Pa}$	$\rho =0.002300\mathrm{kg}{m}^{-3}$, $u=0m{s}^{-1}$, $p=43.772\mathrm{Pa}$

Table 2. Table 2 : Reynolds numbers based on the jet flow velocity computed at t = 23.46 ​ μ ​ s 𝑡 23.46 𝜇 s t=23.46\ \mathrm{\mu s} using pore-resolved numerical simulations.

Aspect ratio of pore geometry	$V_s=279.4m{s}^{-1}$	$V_s=335.9m{s}^{-1}$	$V_s=367.1m{s}^{-1}$
$2/7\times 2/7$	7.3	13.50	17.96
$8\%\times 100\%$	16.43	23.75	27.38

Table 3. Table 3 : Pore-resolved numerical simulations: computational domains and mesh sizes.

Simulation description	Pore geometry	Computational	Mesh size
		domain (mm)
Evaluation of a candidate equivalent pore (Section 4.2)	$2/7\times 2/7$	$4\times 0.5\times 0.5$	$800\times 182\times 182$
Evaluation of a candidate equivalent pore (Section 4.2)	$20\%\times 40\%$	$4\times 0.5\times 0.5$	$800\times 200\times 200$
Evaluation of a candidate equivalent pore (Section 4.2)	$10\%\times 80\%$	$4\times 0.5\times 0.5$	$800\times 200\times 200$
Evaluation of a candidate equivalent pore (Section 4.2)	$8\%\times 100\%$	$4\times 0.5$	$1,600\times 200$
Shock-porous-wall interaction (Section 4.3)	$2/7\times 2/7$	$10\times 0.5\times 0.5$	$2,000\times 84\times 84$
Shock-porous-wall interaction (Section 4.3)	$8\%\times 100\%$	$10\times 0.5$	$2,000\times 200$

Equations

$$\frac{\partial \mathbf{W}}{\partial t}+\nabla \cdot \mathcal{F}(\mathbf{W})=\nabla \cdot \mathcal{G}(\mathbf{W}),$$

$$\bm{W}=\begin{pmatrix}\rho\\[11.38109pt] \rho\bm{v}\\[11.38109pt] E\end{pmatrix},\quad\mathcal{F}(\bm{W})=\begin{pmatrix}\rho\bm{v}\\[5.69054pt] \rho\bm{v}\otimes\bm{v}+p\mathcal{I}\\[7.11317pt] \bigl{(}E+p\bigr{)}\bm{v}\end{pmatrix},\quad\mathrm{and}\quad\mathcal{G}(\bm{W})=\begin{pmatrix}0\\[7.11317pt] {\uptau}\\[4.2679pt] {\uptau}\bm{\cdot}\bm{v}-\bm{q}\end{pmatrix},$$

$$\mathbf{v}=(u,v,w{)}^T\text{and}E=\rho e+\frac{1}{2}\rho (u^2+v^2+w^2),$$

$$\text{uptau}=\mu ({\nabla }^T\mathbf{v}+\nabla \mathbf{v})-\frac{2}{3}\mu (\nabla \cdot \mathbf{v})\mathcal{I}\mathrm{and}\mathbf{q}=-\kappa \nabla T,$$

$$p=\rho RT\mathrm{and}e=\frac{R}{\gamma -1}T,$$

$$\bar{\varphi }=\frac{{\int }_{\Omega }\varphi d\Omega }{|\Omega |},$$

$$\begin{array}{rl}\frac{\partial \bar{\rho }}{\partial t}+\frac{\partial \overline{\rho u}}{\partial x}& =0,\\ \frac{\partial \overline{\rho u}}{\partial t}+\frac{\partial \overline{\rho uu+p}}{\partial x}& =\frac{\partial \overline{{\tau }_{xx}}}{\partial x}+\underset{\text{to be modeled}}{\underbrace{\left[\frac{\overline{\partial {\tau }_{xy}}}{\partial y}+\frac{\overline{\partial {\tau }_{xz}}}{\partial z}\right]}},\\ \frac{\partial \overline{E}}{\partial t}+\frac{\partial \overline{(E+p)u}}{\partial x}& =\frac{\partial \overline{{\tau }_{xx}u-q_x}}{\partial x}+\underset{\text{to be modeled}}{\underbrace{\left[\frac{\overline{\partial {\tau }_{xy}u}}{\partial y}+\frac{\overline{\partial {\tau }_{xz}u}}{\partial z}\right]}},\end{array}$$

$${\tau }_{xx}=\mu \left(\frac{4}{3}\frac{\partial u}{\partial x}-\frac{2}{3}\frac{\partial v}{\partial y}-\frac{2}{3}\frac{\partial w}{\partial z}\right),{\tau }_{xy}=\mu \left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right),{\tau }_{xz}=\mu \left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right).$$

$$q_x=-\kappa \frac{\partial T}{\partial x}.$$

$$\overline{{\tau }_{xx}}=\frac{4}{3}\mu \frac{\partial \overline{u}}{\partial x}\text{ and }\overline{{\tau }_{xx}u-q_x}=\frac{4}{3}\mu \overline{\frac{\partial u}{\partial x}u}+\kappa \frac{\partial \overline{T}}{\partial x},$$

$$\begin{array}{rl}\frac{\partial \bar{\rho }}{\partial t}+\frac{\partial \overline{\rho u}}{\partial x}& =0,\\ \frac{\partial \overline{\rho u}}{\partial t}+\frac{\partial \overline{\rho uu+p}}{\partial x}& =\frac{\partial }{\partial x}\left(\frac{4}{3}\mu \frac{\partial \overline{u}}{\partial x}\right)+\underset{\text{to be modeled}}{\underbrace{\left[\frac{\overline{\partial {\tau }_{xy}}}{\partial y}+\frac{\overline{\partial {\tau }_{xz}}}{\partial z}\right]}},\\ \frac{\partial \overline{E}}{\partial t}+\frac{\partial \overline{(E+p)u}}{\partial x}& =\frac{\partial }{\partial x}\left(\frac{4}{3}\mu \overline{\frac{\partial u}{\partial x}u}+\kappa \frac{\partial \overline{T}}{\partial x}\right)+\underset{\text{to be modeled}}{\underbrace{\left[\frac{\overline{\partial {\tau }_{xy}u}}{\partial y}+\frac{\overline{\partial {\tau }_{xz}u}}{\partial z}\right]}}.\end{array}$$

$$\alpha =\frac{|\Omega {|}^{\text{pore}}}{|\Omega |}$$

$${\mathbf{F}}_{ij}=-{\mathbf{F}}_{ji}=|{\partial \mathcal{C}}_{ij}|{\mathbf{F}}^{\text{Roe}}({\mathbf{W}}_i,{\mathbf{W}}_j,n_{ij}),$$

$${\mathbf{F}}_{ij}=|\partial {\mathcal{C}}_{ij}|{\mathbf{F}}^{\text{imp-wall}}({\mathbf{W}}_i,n_{ij})\text{ and }{\mathbf{F}}_{ji}=|\partial {\mathcal{C}}_{ji}|{\mathbf{F}}^{\text{imp-wall}}({\mathbf{W}}_j,n_{ji}),$$

$$\begin{array}{rl}{\mathbf{F}}_{ij}^{\text{ave}}& =|\partial {\mathcal{C}}_{ij}|\left(\alpha {\mathbf{F}}^{\text{Roe}}({\overline{\mathbf{W}}}_i,{\overline{\mathbf{W}}}_j,n_{ij})+(1-\alpha ){\mathbf{F}}^{\text{imp-wall}}({\overline{\mathbf{W}}}_i,n_{ij})\right),\\ {\mathbf{F}}_{ji}^{\text{ave}}& =|\partial {\mathcal{C}}_{ji}|\left(\alpha {\mathbf{F}}^{\text{Roe}}({\overline{\mathbf{W}}}_j,{\overline{\mathbf{W}}}_i,n_{ji})+(1-\alpha ){\mathbf{F}}^{\text{imp-wall}}({\overline{\mathbf{W}}}_j,n_{ji})\right),\end{array}$$

$$\begin{array}{r}{\mathbf{G}}_{ij}^{\text{ave}}=|\partial {\mathcal{C}}_{ij}|\left(\alpha \underset{{\mathbf{G}}_{ij}^{\text{pore}}}{\underbrace{\left(\begin{array}{c}0\\ \frac{4}{3}\mu \frac{{\overline{u}}_j-{\overline{u}}_i}{{\delta }_x}\\ \frac{4}{3}\mu \left(\frac{{\overline{u}}_j-{\overline{u}}_i}{{\delta }_x}\right)\left(\frac{{\overline{u}}_i+{\overline{u}}_j}{2}\right)+\kappa \frac{{\overline{T}}_j-{\overline{T}}_i}{{\delta }_x}\end{array}\right)}}+(1-\alpha )\underset{{\mathbf{G}}_{ij}^{\text{imp-wall}}}{\underbrace{\left(\begin{array}{c}0\\ \frac{4}{3}\mu \frac{{\overline{u}}_j^{\text{be}}-{\overline{u}}_i}{{\delta }_x}\\ \frac{4}{3}\mu \left(\frac{{\overline{u}}_j^{\text{be}}-{\overline{u}}_i}{{\delta }_x}\right)\left(\frac{{\overline{u}}_i+{\overline{u}}_j^{\text{be}}}{2}\right)+\kappa \frac{{\overline{T}}_j^{\text{be}}-{\overline{T}}_i}{{\delta }_x}\end{array}\right)}}\right),\\ {\mathbf{G}}_{ji}^{\text{ave}}=|{\partial \mathcal{C}}_{ij}|\left(\alpha \underset{{\mathbf{G}}_{ji}^{\text{pore}}}{\underbrace{\left(\begin{array}{c}0\\ \frac{4}{3}\mu \frac{{\overline{u}}_i-{\overline{u}}_j}{{\delta }_x}\\ \frac{4}{3}\mu \left(\frac{{\overline{u}}_i-{\overline{u}}_j}{{\delta }_x}\right)\left(\frac{{\overline{u}}_i+{\overline{u}}_j}{2}\right)+\kappa \frac{{\overline{T}}_i-{\overline{T}}_j}{{\delta }_x}\end{array}\right)}}+(1-\alpha )\underset{{\mathbf{G}}_{ji}^{\text{imp-wall}}}{\underbrace{\left(\begin{array}{c}0\\ \frac{4}{3}\mu \frac{{\overline{u}}_i^{\text{be}}-{\overline{u}}_j}{{\delta }_x}\\ \frac{4}{3}\mu \left(\frac{{\overline{u}}_i^{\text{be}}-{\overline{u}}_j}{{\delta }_x}\right)\left(\frac{{\overline{u}}_i^{\text{be}}+{\overline{u}}_j}{2}\right)+\kappa \frac{{\overline{T}}_i^{\text{be}}-{\overline{T}}_j}{{\delta }_x}\end{array}\right)}}\right),\end{array}$$

$$u(y,z)=C\left(r^2-y^2-z^2\right),$$

$$\overline{u}=\frac{{\int }_{\Omega }udydz}{|\Omega |}=\frac{{\int }_{{\Omega }^{\text{pore}}}udydz}{|\Omega |}=C\frac{|{\Omega }^{\text{pore}}|}{|\Omega |}\frac{r^2}{2}=C\alpha \frac{r^2}{2}.$$

$$\begin{array}{rl}\overline{\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial z^2}}& =-8\frac{\overline{u}}{r^2},\\ \overline{u\left(\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial z^2}\right)}& =-8\frac{{\overline{u}}^2}{\alpha r^2},\\ \overline{{\left(\frac{\partial u}{\partial y}\right)}^2+{\left(\frac{\partial u}{\partial z}\right)}^2}& =8\frac{{\overline{u}}^2}{\alpha r^2}.\end{array}$$

$$\begin{array}{rl}\overline{\frac{\partial {\tau }_{xy}}{\partial y}}+\overline{\frac{\partial {\tau }_{xz}}{\partial z}}& =-8\mu \frac{\overline{u}}{r^2},\\ \overline{\frac{\partial {\tau }_{xy}u}{\partial y}}+\overline{\frac{\partial {\tau }_{xz}u}{\partial z}}& =0.\end{array}$$

$$u(y,z)=C\frac{16r^2}{{\pi }^3}\sum _{i=1,3,\mathrm{5...}}^{\infty }\left[(-1{)}^{(i-1)/2}\left(1-\frac{\cosh (i\pi z/2r)}{\cosh (i\pi /2)}\right)\frac{\cos (i\pi y/2r)}{i^3}\right],$$

$$\overline{u}$$

$$\begin{split}\overline{\frac{\partial^{2}u}{\partial y^{2}}+\frac{\partial^{2}u}{\partial z^{2}}}&=-\frac{3}{\zeta}\frac{\overline{u}}{r^{2}}\,,\\ \overline{u\left(\frac{\partial^{2}u}{\partial y^{2}}+\frac{\partial^{2}u}{\partial z^{2}}\right)}+\overline{\left(\frac{\partial u}{\partial y}\Big{)}^{2}+\Big{(}\frac{\partial u}{\partial z}\right)^{2}}&=0\,,\end{split}$$

$$\zeta =1-\frac{192}{{\pi }^5}\sum _{i=1,3,\mathrm{5...}}^{\infty }\frac{\tanh (i\pi /2)}{i^5}\approx \mathrm{0.421731044865.}$$

$$\begin{array}{rl}\overline{\frac{\partial {\tau }_{xy}}{\partial y}}+\overline{\frac{\partial {\tau }_{xz}}{\partial z}}& =-\frac{3}{\zeta }\mu \frac{\overline{u}}{r^2},\\ \overline{\frac{\partial {\tau }_{xy}u}{\partial y}}+\overline{\frac{\partial {\tau }_{xz}u}{\partial z}}& =0.\end{array}$$

$$u(y,z)=C\left(r^2-y^2\right),$$

$$\overline{u}=\frac{{\int }_{\Omega }udydz}{|\Omega |}=\frac{{\int }_{{\Omega }^{\text{pore}}}udydz}{|\Omega |}=C\frac{|{\Omega }^{\text{pore}}|}{|\Omega |}\frac{2r^2}{3}=C\alpha \frac{2r^2}{3}.$$

$$\begin{array}{rl}\overline{\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial z^2}}& =-3\frac{\overline{u}}{r^2},\\ \overline{u\left(\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial z^2}\right)}& =-3\frac{{\overline{u}}^2}{\alpha r^2},\\ \overline{{\left(\frac{\partial u}{\partial y}\right)}^2+{\left(\frac{\partial u}{\partial z}\right)}^2}& =3\frac{{\overline{u}}^2}{\alpha r^2}.\end{array}$$

$$\begin{array}{rl}\overline{\frac{\partial {\tau }_{xy}}{\partial y}}+\overline{\frac{\partial {\tau }_{xz}}{\partial z}}& =-3\mu \frac{\overline{u}}{r^2},\\ \overline{\frac{\partial {\tau }_{xy}u}{\partial y}}+\overline{\frac{\partial {\tau }_{xz}u}{\partial z}}& =0.\end{array}$$

$$D(x)=\exp \left[-\frac{1}{2}{\left(\frac{x-x_0}{\sigma }\right)}^2\right],$$

Full text

1 Introduction

Porous walls and membranes appear in a wide range of scientific and engineering applications. These include, to name only a few, filtration processes [27, 11, 21] to separate particulates or molecules from a bulk fluid, the cooling of gas turbines via multiperforated plates [1, 24, 23, 25], light weight canopies for low density supersonic decelerators [18, 17, 9, 15, 14], and windbreaks or shelterbelts [16, 29, 37] to reduce the relative flow speed and maintain stability. For most of these large-scale engineering applications, pore-resolved CFD computations remain unpractical, if not unaffordable. Therefore, the high-precision and yet computationally efficient treatment of porous wall boundary conditions in CFD computations is crucial for enabling flow simulations that further the understanding of the effects of porous materials and can assist in the development of membrane-related processes and technologies.

The modeling of fluid flow through a porous wall or membrane has been inspired by a variety of investigations of flows in porous media. In [18, 36, 34, 9, 35], a modified Darcy’s law was coupled with Ergun’s equation to construct an internal pressure jump condition for the incompressible or low-Mach Navier-Stokes equations and account for the material porosity of a canopy in the simulation of parachute fluid-structure interactions. In [40] and [37], the flow associated with a windbreak similar to a plant canopy was modeled by adding a local, momentum extraction, source term near the windbreak. The pressure-loss coefficient in the source term was calibrated using the pressure drop across the windbreak, as well as the flow momentum measured in related wind tunnel experiments [16]. In [27] and [11], the analysis of the cross-flow filtration was performed by numerically simulating the transport of solutes or particles in a tube with porous walls using the Navier-Stokes equations and applying the slip wall boundary condition, as advocated in [3]. Notably, all of the aforementioned modeling approaches are limited to the incompressible flow regime.

Only a few approaches for modeling the effects of permeability on high Mach number compressible flows have been documented so far. Those approaches have combined Darcy’s law with the continuity equation [26], or derived compressible variants of Darcy’s law from the continuity equation and the force balance equation between drag and the pressure drop across an element of a porous material [32, 31]. For flow computations associated with film-cooling in gas turbines operating at high-speed, the following works are noteworthy. In [1], the treatment of the porous wall boundary condition was performed using a mass inflow condition with white noise in the mass flow rate. In [23], a homogeneous model was proposed for a multiperforated plate where the homogenized mass flux is related to a constant discharge coefficient estimated from pore-resolved numerical simulations.

Here, an alternative homogenization approach is proposed for the treatment of porous wall boundary conditions in the numerical simulation of compressible viscous flows through permeable surfaces. Unlike the aforementioned approaches, it does not require prescribing a mass flow rate nor introduces an empirical discharge coefficient in the computational model. Instead, it models the inviscid flux through a porous wall as a weighted average of the inviscid flux at an impermeable surface and that through pores. To account for the loss of friction along the pore boundaries, it introduces a body force term in the governing equations that depends on the thickness of the porous wall and the concept of an “equivalent pore” – that is, a single pore whose geometry and dimensions are such that when inserted at the center of a unit cell, the void fraction of the cell matches that of the fabric of interest. The feasibility of such a concept is demonstrated using low-speed permeability experimental data for the fabric of the Mars Science Laboratory parachute canopy. The overall homogenization approach can be implemented in any finite difference, finite volume, or finite element semi-discretization scheme. It is illustrated with a series of supersonic porous flow problems with Martian atmosphere free-stream conditions and verified using pore-resolved numerical simulations.

The remainder of the paper is organized as follows. First, the governing equations are specified in Section 2. Next, the proposed homogenization approach is presented in Section 3 and its verification is performed in Section 4. Finally, conclusions are offered in Section 5.

2 Governing equations

Throughout this paper, flows are assumed to be compressible as well as viscous and are governed by the Navier-Stokes equations. These can be written in conservation form as follows

[TABLE]

where $\bm{W}$ denotes the vector of the conservative variables describing the fluid state, $t$ denotes time, and $\mathcal{F}(\bm{W})$ and $\mathcal{G}(\bm{W})$ are the inviscid and viscous flux tensor functions of the conservative variables, respectively. Specifically,

[TABLE]

where $\rho$, $\bm{v}$, and $E$ denote the density, velocity, and total energy per unit volume of the fluid, respectively, and $\otimes$ denotes the Hadamard product. The fluid velocity and total energy per unit volume can be written as

[TABLE]

where the superscript $T$ designates the transpose and $e$ denotes the specific (i.e., per unit of mass) internal energy. In the physical inviscid flux tensor $\mathcal{F}$, $p$ is the static pressure and $\mathcal{I}\in\mathbb{R}^{3}$ is the identity matrix. In the viscous flux tensor $\mathcal{G}$, ${\uptau}$ and $\bm{q}$ denote the viscous stress tensor and heat flux vector, respectively, and are defined as follows

[TABLE]

where $\mu$ is the dynamic shear viscosity of the fluid, $\kappa$ is its thermal conductivity, and $T$ is its temperature.

For simplicity, but without any loss of generality, the system of equations (1) is closed by assuming that the gas is ideal, calorically perfect, and therefore governed by

[TABLE]

where $R$ and $\gamma$ are the gas constant and specific heat ratio, respectively.

3 Homogenized flux-body force approach

Consider as a backdrop for the homogenization approach described herein the three-dimensional, porous flow model problem graphically depicted in Fig. 1. The computational fluid domain is a three-dimensional, orthogonal box with a square cross section (the generalization to an arbitrarily-shaped computational fluid domain is straightforward). It contains at its center a porous wall of a finite thickness (shown in blue color) whose edges are flushed with four of the six surface boundaries of the computational fluid domain. The free-stream velocity of the flow is perpendicular to the surface of the wall and periodic boundary conditions are prescribed in both transverse directions.

Assuming that the porous wall is reasonably thin, it can be modeled as a surface denoted here by $\Omega$. This surface is discretized by a mesh whose typical element size is much larger than the pore size. Hence, given a scalar or vector quantity $\phi$ defined in $\Omega$, a homogenized semi-discrete counterpart $\bar{\phi}$, where homogenization is performed in both transverse directions and designated by an overline, is defined here as

[TABLE]

where $|\Omega|$ denotes the area of the surface $\Omega$.

Due to symmetry (see Fig. 1), both homogenized transverse velocity components $\bar{v}$ and $\bar{w}$ are zero. For this reason, the homogenized version of Eq. 1 inside the porous wall can be written as

[TABLE]

where

$\tau_{xx}$, $\tau_{xy}$, and $\tau_{xz}$ are the components of the viscous stress tensor ${\uptau}$ given by

[TABLE]

• 2.

$q_{x}$ is the $x$-component of the heat flux vector $\bm{q}$ whose expression is given by

[TABLE]

• 3.

The averaged terms $\overline{\partial\tau_{xy}/\partial y}+\overline{\partial\tau_{xz}/\partial z}$ and $\overline{\partial\tau_{xy}u/\partial y}+\overline{\partial\tau_{xz}u/\partial z}$ are shown between brackets to emphasize that: while they vanish in the cross sections away from the porous wall – due to the fact that they contain only partial derivatives with respect to the directions of homogenization $y$ and $z$ – they represent the averaged friction loads on the pore boundary and therefore require subgrid scale modeling; and hence they will be reintroduced in Section 3.3 in the form of a body force model.

From (8), (9), and the periodic nature of the prescribed boundary conditions (see Fig. 1), it follows that the expressions for the averaged, wall normal, viscous flux terms $\partial\overline{\tau_{xx}}/\partial x$ and $\partial\left(\overline{\tau_{xx}u-q_{x}}\right)/\partial x$ can be simplified as follows

[TABLE]

which transforms Eq. 7 into the following, more specific, homogenized governing equations inside the porous wall

[TABLE]

For the sake of clarity, the homogenization approach proposed herein is described in the next two subsections not only in the backdrop of the model porous flow problem illustrated in Fig. 1 but also in the context of a cell-centered Finite Volume (FV) approximation method, where the gradients underlying the viscous fluxes are semi-discretized by a finite difference scheme and for simplicity, but without any loss of generality, the CFD mesh is assumed to be structured and body-fitted. In this case, each computational cell ${\mathcal{C}}_{i}$ is a primal element of the CFD mesh, $i$ is the center point of ${\mathcal{C}}_{i}$ where the semi-discrete fluid state vector is defined, and the boundary surface (cell facet) $\partial{\mathcal{C}}_{ij}$ shared by two cells $\mathcal{C}_{i}$ and $\mathcal{C}_{j}$, where $j$ is such that the line $(i,j)$ traverses the porous wall, lies along the surface of this wall (see Fig. 2). However, as will be noted wherever appropriate within the next two subsections, the generalizations of this approach to a vertex-based FV approximation method where the viscous fluxes are semi-discretized by a Finite Element (FE) method and to unstructured meshes as well as an arbitrary orientation of the porous wall with respect to the global reference frame, are straightforward. Furthermore, the homogenization approach described below is not limited to body-fitted CFD methods: it is equally applicable to embedded/immersed boundary methods. In particular, it is well-suited for the FInite Volume method with Exact two-material Riemann problems (FIVER) [38, 19, 22, 12, 4] (see also Appendix A). As a matter of fact, all numerical results with the porosity model reported in Section 4 are obtained using non body-fitted CFD meshes and the Embedded Boundary Method (EBM) FIVER. In short, the proposed homogenization approach for the treatment of porous wall boundary conditions can be incorporated in any semi-discretization or discretization method.

3.1 Homogenized inviscid flux

Let $\alpha$ denote the void fraction or porosity of the porous wall. It can be expressed as

[TABLE]

where $|\Omega^{\textrm{pore}}|$ denotes the total area occupied by the pores in the surface $\Omega$.

Consider first the two following extreme cases:

$\alpha=1$, which corresponds to the no-wall configuration. In this case, the inviscid fluxes through the (nonexistent) porous wall can be approximated as usual, for example, by an approximate Riemann solver such as Roe’s solver [30]. This can be written as (see Fig. 2)

[TABLE]

where $\bm{F}_{ij}$ is the numerical inviscid flux function associated with the inviscid flux vector in the $x$-direction, $\bm{F}$, $n_{ij}$ is the unit outward normal to $\partial{\mathcal{C}}_{ij}$ – in this case, $n_{ij}$ is parallel to the $x$-direction – and $|\partial{\mathcal{C}}_{ij}|$ denotes the area of $\partial{\mathcal{C}}_{ij}$.

• 2.

$\alpha=0$, which corresponds to the impermeable wall configuration. In this case, the inviscid fluxes can also be approximated as usual – that is, as wall boundary fluxes as follows (see Fig. 2)

[TABLE]

where $|\partial{\mathcal{C}}_{ji}|=|\partial{\mathcal{C}}_{ij}|$, $n_{ji}=-n_{ij}$, and $\bm{F}^{\textrm{{imp-wall}}}(\bm{W}_{i},n_{ij})$ and $\bm{F}^{\textrm{{imp-wall}}}(\bm{W}_{j},n_{ji})$ are computed using the physical inviscid flux tensor $\mathcal{F}$, $n_{ij}$, and the standard, impermeable wall boundary conditions.

Hence, the main idea here is to compute, for any nontrivial porosity value $0<\alpha<1$, the homogenized, numerical inviscid flux functions through the porous wall $\bm{F}^{\textrm{ave}}_{ij}$ and $\bm{F}^{\textrm{ave}}_{ji}$ as convex combinations of (13) and (14) using however averaged values of the fluid state vectors – that is,

[TABLE]

where averaging is performed as in (6). The second-order extension of the above approximations using the classical Monotonic Upstream-centered Scheme for Conservation Laws (MUSCL) is straightforward.

The above computation (15) is illustrated in Fig. 2, where the computational cells adjacent to the porous wall are divided into two parts: one part attached to a hole (colored in blue) representing the pores, collectively; and another part attached to the wall (colored in red). Accordingly, the flux through the porous wall $\bm{F}^{\textrm{ave}}_{ij}$ consists of the flux through the hole (13) and that at the impermeable part of the wall (see Eq. 14).

Note that by construction, the proposed homogenization of the numerical inviscid flux function (15) associated with the inviscid flux vector in the $x$-direction $\bm{F}$ conserves mass and energy; it also guarantees that the impermeable part of the numerical flux has zero mass and energy components.

The generalization of the inviscid flux homogenization procedure described above to the case of a vertex-based FV method with control volumes (dual cells) ${\mathcal{C}}_{i}$ – in which case a wall boundary grid point $i$ and not a wall boundary cell facet $\partial{\mathcal{C}}_{ij}$ lies on the porous wall – is simply accomplished using the concept of half dual cells. The generalization of this procedure to an arbitrary direction of the outward normal to the wall $n_{ij}$ and to unstructured meshes is straightforward, because most if not all approximate Riemann solvers are essentially one-dimensional local solvers that propagate information in the direction normal to the considered facet of the control volume boundary.

3.2 Homogenized viscous flux

The spirit of the approach described in Section 3.1 for the numerical inviscid fluxes is followed here to model each numerical viscous flux in the wall normal direction (which in this case is the $x$-direction), $\bm{G}_{ij}^{\textrm{ave}}$ and $\bm{G}_{ji}^{\textrm{ave}}$, as a convex combination of the numerical viscous flux at an impermeable surface and that through pores computed using averaged fluid state vectors. Recalling that the terms shown between brackets in Eq. 11 vanish because they contain derivatives in the homogenized directions, choosing to approximate the gradients in the viscous fluxes by, for example, the second-order central finite differencing, and assuming that the dynamic shear viscosity $\mu$ and fluid thermal conductivity $\kappa$ are constant, the homogenization of the numerical viscous flux functions through the porous wall can be carried out as follows (see Fig. 2)

[TABLE]

where $\delta_{x}$ denotes the mesh spacing in the $x$-direction and $\overline{u}^{\,\textrm{be}}_{i}$, $\overline{u}^{\,\textrm{be}}_{j}$, $\overline{T}^{\,\textrm{be}}_{i}$, and $\overline{T}^{\,\textrm{be}}_{j}$ are nodal velocity and temperature values populated using any preferred extrapolation approaches that enforce the no-slip and adiabatic wall boundary conditions in cell-centered FV approximation methods.

As stated earlier, the vanishing terms $\overline{\partial\tau_{xy}/\partial y}+\overline{\partial\tau_{xz}/\partial z}$ and $\overline{\partial\tau_{xy}u/\partial y}+\overline{\partial\tau_{xz}u/\partial z}$ of Eq. 11 must be modeled as they represent the averaged friction loads (per unit volume) at a pore boundary (or at an expanded jet stream boundary). The evaluation of these terms requires the availability of a profile of the component of the flow velocity vector normal to the porous wall – in this case, $u(y,z)$ – that is fully resolved inside a pore: hence, it requires subgrid scale modeling. Such modeling is inspired here by the Hagen-Poiseuille flow [33], which in this case is homogeneous in the $x$-direction. It starts by postulating: an equivalent pore – that is, a single pore whose shape and dimensions are such that when inserted at the center of a unit cell, the void fraction of the cell matches that of the fabric of interest; and a parabolic shape for the profile of the wall normal velocity component in the equivalent pore that depends on its geometry. The subgrid scale modeling ends by deriving an analytical expression of the form of a body force term for the homogenized equations of dynamic equilibrium (11). Then, the body force term, which represents the friction loss along pore boundaries, is approximated by any preferred numerical procedure.

To this end, circular, square, and gapped shapes pores (see Fig. 3) are examined next. For each of these geometries, which can be considered as an idealization of a family of more complex geometries, a velocity profile $u(y,z)$ inside the pore is postulated and the form of the corresponding body force term is derived.

At this point, it is noted that the extension of the homogenization of the numerical viscous flux functions (16) to the case of a vertex-based FV approximation method, where the viscous fluxes are semi-discretized by a FE method, is straightforward. In such a context, the computational cell ${\mathcal{C}}_{i}$ becomes the dual cell associated with the grid point $i$ of the CFD mesh, which can be structured or arbitrarily unstructured. If this mesh is body-fitted, all wall boundary grid points lie on the wall boundary. In this subcase, $|{\mathcal{C}}_{i}|$ is replaced by the volume of a primal element $B^{e}$ of the CFD mesh that is connected to the grid point $i$; the computation of the term ${\bm{G}_{ij}^{\textrm{pore}}}$ is unchanged except that the constant gradients in this term are computed by differentiating the linear shape functions of the element $B^{e}$; and summation is performed over all elements $B^{e}$ connected to the grid point $i$. As for the term ${\bm{G}_{ij}^{\textrm{imp-wall}}}$, it simplifies to the standard computation of the wall boundary viscous fluxes – that is, no extrapolation is needed; and similarly, summation is performed over all elements $B^{e}$ connected to the grid point $i$. On the other hand, if the CFD mesh is not body-fitted, a preferred embedded/immersed boundary method can be used as is or tailored to compute both terms ${\bm{G}_{ij}^{\textrm{pore}}}$ and ${\bm{G}_{ij}^{\textrm{imp-wall}}}$ (for example, see Appendix A), whether the FV approximation method is cell-centered or vertex-based.

3.2.1 Circular equivalent pore

For a circular equivalent pore such as that shown in Fig. 3(a), the Hagen-Poiseuille flow theory leads to the following velocity profile inside the pore [39]

[TABLE]

where $r$ is the radius of the circular geometry and $C$ is a free parameter that will be shown not play an explicit role in the computation of the body force term. From (6), (12), and (17), it follows that in this case, the averaged flow velocity is

[TABLE]

From (17) and (18), it follows that

[TABLE]

Substituting the results (19) into the expressions (8) leads to the following subgrid scale model for the friction loss along circular pore boundaries

[TABLE]

3.2.2 Square equivalent pore

For a square equivalent pore such as that shown in Fig. 3(b), the Hagen-Poiseuille flow theory leads to the following velocity profile inside the pore [39, p. 113]

[TABLE]

where $2r$ is the side length of the square geometry and $C$ is a free parameter that will be shown not play an explicit role in the computation of the body force term. From (6), (12), and (21), it follows that in this case, the averaged flow velocity is given by

[TABLE]

From (21) and (22), it follows that

[TABLE]

where

[TABLE]

Substituting (23) into (8) leads to the following subgrid scale model for the friction loss along the boundaries of a square pore

[TABLE]

3.2.3 Gapped shape equivalent pore

For the gapped shape equivalent pore shown in Fig. 3(c), the Hagen-Poiseuille flow theory leads to the following velocity profile inside the pore [39]

[TABLE]

where $2r$ is in this case the width of the gap and $C$ is a free parameter that will be shown not to play an explicit role in the computation of the body force term. From (6), (12), and (25, it follows that in this case, the averaged flow velocity can be written as

[TABLE]

From (25) and (26), it follows that

[TABLE]

Substituting the results (27) into (8) yields the following subgrid scale model for the friction loss along the boundaries of a gapped shape pore

[TABLE]

3.3 Body force term

From the second of each of Eqs. (20), (24), and (28), it follows that whatever symmetric pore geometry is assumed, the resulting friction loss contributes only to a momentum loss. In particular, it does not contribute to a total energy loss, because friction converts the lost kinetic energy into internal energy. For $\alpha=8\%$, the momentum loss for a circular pore with $\displaystyle r=\sqrt{{\alpha|\Omega|}/{\pi}}$ represented by the form of the body force term derived in (20) is about $\displaystyle-314.16\mu(\overline{u}/|\Omega|)$, that for a square pore with $\displaystyle r=\sqrt{{\alpha|\Omega|}/{4}}$ represented by the counterpart form (24) is roughly $\displaystyle-355.68\mu(\overline{u}/|\Omega|)$, and that due to a gapped shape pore with $\displaystyle r={(\alpha\sqrt{|\Omega|})}/{2}$ and the form of its body force term (28) is around $\displaystyle-1875\mu(\overline{u}/|\Omega|)$. This observation suggests that the smaller is the aspect ratio of a pore, the smaller is the associated momentum loss or dissipation.

In any case, the derived form of the body force term is multiplied by the following Gaussian approximation of the Dirac function in order to locate it along the pore boundaries

[TABLE]

where $x_{0}$ is the position of the porous wall – in this case along the $x$-direction – $\sigma=\eta_{f}/(2\pi)$, and $\eta$ is the thickness of the porous wall (note that the integral of $D(x)$ along the thickness of a porous wall is equal to the thickness of the wall). Then, the following body force (source) term is added to the homogenized equations of dynamic equilibrium (11)

[TABLE]

where

$\eta_{f}$ is a thickness correction factor that transforms the thickness of the porous wall into an “effective” thickness in order to account for the expanded jet stream through a pore, which also contributes to the momentum loss.

• 2.

The term $(\overline{\partial\tau_{xy}/\partial y}+\overline{\partial\tau_{xz}/\partial z})$ is given in the first of Eqs. (20) in the case of a circular pore, the first of Eqs. (24) in the case of a square pore, and the first of Eqs. (28) in the case of a gapped-shape pore.

• 3.

In each of the first of Eqs. (20), the first of Eqs. (24), and the first of Eqs. (28), $\bar{u}$ is not computed using (18), (22), and (26), respectively. Instead, $\bar{u}$ is explicited computed using the averaging definition (6). Consequently, the free parameter $C$ in (17), (21), and (25) is never used explicitly in the computation of the body force term (30).

4 Demonstration and Verification

As stated above, the main interest here is in compressible and particularly supersonic flows through permeable surfaces. A notable application is the prediction of parachute inflation dynamics during Mars landing operations [15, 4, 13]. For this application, a relevant fabric is the PIA-C-7020 Type I fabric with Ripstop weave selected for the parachute canopy of the Mars Science Laboratory. Based on the X-ray microtomography results published in [28], the porosity of this fabric (see Fig. 4) is estimated at roughly 8%. Its thickness is around $80\ \mathrm{\mu m}$ [20]. In both in-plane directions, the period of the weave pattern is about $500\ \mathrm{\mu m}$.

To the best of knowledge of the authors, only low-speed permeability test data is available for the PIA-C-7020 Type I fabric with Ripstop weave (and for that matter, any other fabric material). Hence, the proposed homogenization approach for the treatment of porous wall boundary conditions is demonstrated and verified here for supersonic flows as follows:

First, a most-suitable equivalent pore for the complex Ripstop weave pattern is determined by exploring the potential of various candidate pore shapes. For each such candidate: a series of pore-resolved, low-speed, unsteady numerical simulations is performed; and the computed variation of the averaged velocity through the fabric with the difference between the prescribed upstream and downstream pressure values, which defines the permeability of the fabric, is compared to its counterpart obtained from available low-speed permeability test data. Among all explored candidate pore geometries, the most-suitable one is identified as that which leads to the closest value of the experimentally measured permeability of the fabric.

• 2.

Next, supersonic, pore-resolved, unsteady numerical simulations are performed for the identified most-appropriate pore geometry and for another test geometry in order to generate high-speed reference flow results in lieu of high-speed permeability test data. Counterpart unsteady simulations using the proposed homogenization approach are also carried out. The obtained numerical results are compared to their reference values.

4.1 General problem setup

Given the periodicity of the weave pattern, all aforementioned numerical simulations are performed at the coupon level. Specifically, all unsteady simulations are carried out for a small square piece of fabric of edge size equal to $500\ \mathrm{\mu m}$, which corresponds to the period of the weave pattern. For simplicity, all flexibility effects of the coupon are neglected. For the pore-resolved numerical simulations, the weave pattern is simplified by lumping all gaps between spun fibers into a single square or rectangular hole that is placed at the center of the coupon and sized so that the void fraction $\alpha$ of the material is preserved.

The explored equivalent pores are shown in Fig. 5. All have the same porosity $\alpha=8\%$, but different aspect ratios – specifically, $2/7\times 2/7$, $20\%\times 40\%$, $10\%\times 80\%$, and $8\%\times 100\%$. The fabric coupon is placed at the center of the computational fluid domain, as in Fig. 1. Periodic boundary conditions are prescribed in both transverse directions, but various initial and uniform upstream/downstream boundary conditions are considered in order to mimic the experimental setup described next.

4.2 Exploration of candidate equivalent pores

The permeability of the PIA-C-7020 Type I fabric with Ripstop weave was measured and reported in [5]. A series of tests were conducted using a Textest Instruments FX 3300 Labotester III Air Permeability Tester, which consists of a test chamber with a vacuum pump. In all tests, the fabric sample was clamped over the test head opening and the upstream flow pressure was set to $p=94,625\ \mathrm{Pa}$. The temperature of the laboratory where the experiment was conducted was measured to be $T=296.5\ \mathrm{K}$ and the following air properties were assumed

[TABLE]

Each different test was performed using still air and a different downstream pressure value that was maintained constant in time using the vacuum pump in the chamber. Then, the permeability of the fabric material was deduced from the experimentally obtained variation of the averaged velocity of the fluid through the fabric with the prescribed difference between the upstream and downstream pressure values.

For each candidate equivalent pore shown in Fig. 5, seven pore-resolved, unsteady numerical simulations mimicking the experimental setup described above are performed using direct numerical simulation (DNS) flow solver Hydrodynamics Adaptive Mesh Refinement Simulator (HAMeRS) [41] as described in Appendix B. In all simulations, the computational fluid domain is chosen to be the $4\ \mathrm{mm}\times 0.5\ \mathrm{mm}\times 0.5\ \mathrm{mm}$ box (dimensions in the $x$-, $y$-, and $z$-directions, respectively) shown in Fig. 1. The upstream (inlet) boundary conditions as well as the initial conditions in the half domain delimited in the $x$-direction by the inlet boundary and the porous wall are set to the uniform flow field defined by $\rho^{\textrm{up}}=1.111\ \mathrm{kg/m}^{3}$, ${\bm{v}}^{\textrm{up}}=0\ \mathrm{m\ s^{-1}}$, and $p^{\textrm{up}}=94,625\ \mathrm{Pa}$. As for the downstream (outlet) boundary conditions and the initial conditions in the other half of the computational fluid domain, they are set in each case to the uniform flow field defined by $\rho^{\textrm{dn}}=1.112\ \mathrm{kg/m}^{3}$ ($1.109\ \mathrm{kg/m}^{3}$, $1.105\ \mathrm{kg/m}^{3}$, $1.099\ \mathrm{kg/m}^{3}$, $1.088\ \mathrm{kg/m}^{3}$, $1.076\ \mathrm{kg/m}^{3}$, and $1.064\mathrm{kg/m}^{3}$), ${\bm{v}}^{\textrm{dn}}=0\ \mathrm{m\ s^{-1}}$, and $p^{\textrm{dn}}=94,525\ \mathrm{Pa}$ ($94,425\ \mathrm{Pa}$, $94,125\ \mathrm{Pa}$, $93,625\ \mathrm{Pa}$, $92,625\ \mathrm{Pa}$, $91,625\ \mathrm{Pa}$, and $90,625\mathrm{Pa}$). Note that these initial conditions are consistent with the temperature $T=296.5\ \mathrm{K}$ of the quiescent air in the laboratory.

At the beginning of each pore-resolved, unsteady numerical simulation characterized by a different uniform outlet pressure field, pressure waves propagate in both directions away from the porous fabric coupon. For each explored pore geometry, the pressure difference across the porous wall instantly drops, the fluid flow near the porous wall accelerates, and a weak jet stream forms downstream. In Fig. 6, the reader can observe that the computed Mach number profiles substantially differ across the considered pore geometries. As expected, the holes with the smaller aspect ratios ($2/7\times 2/7\textrm{ and }20\%\times 40\%$) lead to larger maximum Mach number values (around 0.2), which is consistent with the observation made in Section 3.3 that equivalent pores with smaller aspect ratios lead to smaller momentum dissipation.

In each of the 28 performed pore-resolved numerical simulations, the pressure difference across the porous wall rapidly tends to become quasi-steady and the velocity field on each side of the porous wall converges to a steady-state value around $t=15.8\ \mathrm{\mu s}$. These results are reported in Fig. 7 in the form of velocity-pressure difference curves together with their counterpart experimental values measured in [5] within a confidence interval of $99.7\%$. Here, velocity refers to the averaged quasi-steady velocity through the fabric and pressure difference refers to the averaged quasi-steady pressure difference between both sides of the porous wall.

The permeability velocity through the porous membrane $u^{\textrm{pe}}$ is defined as the ratio of the averaged mass flux across $\Omega$ and the initial upstream density $\rho^{\textrm{up}}$ – that is,

[TABLE]

Both experimental and numerical results reported above are close to the incompressible flow regime, where the density fluctuation can be neglected. Indeed, all performed pore-resolved numerical simulations reveal a density fluctuation of the order of $2\%$. In this case, the definition (32) simplifies to $u^{\textrm{pe}}=\bar{u}$. From this result and Fig. 7, it follows that: each explored pore geometry leads to a numerical prediction of the permeability velocity that matches reasonably well its experimentally measured counterpart; but the geometry labeled as (d) in Fig. 5, which is characterized by the highest aspect ratio ($8\%\times 100\%$ or 12.5:1) 111The geometry labeled as (c) in Fig. 5 delivers similar results because its aspect ratio is very close to that of the geometry labeled as (d)., delivers the closest permeability velocity value to the experimentally measured one and therefore is the most suitable equivalent pore geometry.

It is worth noting that according to [6] and for a given wall permeability, the total pressure drop across a porous wall $(\Delta p)^{\textrm{tot}}$ can be split into two components due to viscous and inertial effects – that is,

[TABLE]

The viscous effect is dominant at low Reynolds number and the inertial effect prevails at high Reynolds number. The study carried out herein falls into the former category ($10\leq Re\leq 400$), where one contribution of the viscous effect is the friction loss due to the shear stress near the surface of each spun fiber. Pore geometries characterized by smaller aspect ratios (e.g., $2/7\times 2/7$ and $20\%\times 40\%$) are also characterized by smaller pore surfaces (e.g., see Fig. 5): hence, they lead to smaller values of $\Delta p^{\textrm{vis}}$ and therefore smaller values of $\Delta p^{\textrm{tot}}$. This explains why the pore geometries with the larger aspect ratios (e.g., $10\%\times 80\%$ and $8\%\times 100\%$) are found in Fig. 7 to lead to numerical predictions of the permeability velocity that match better the experimental data.

4.3 Verification for a supersonic application

Here, the homogenization approach presented in this paper is verified for a series of supersonic instances of the porous flow model problem graphically depicted in Fig. 1. In each instance:

The computational fluid domain is chosen to be the $10\ \mathrm{mm}\times 0.5\ \mathrm{mm}\times 0.5\ \mathrm{mm}$ box (dimensions in the $x$-, $y$-, and $z$-directions, respectively).

• 2.

This domain is divided in two non overlapping subdomains labeled as “Upstream” and “Downstream”: the interface $\Gamma$ between them is located at $250\ \mathrm{\mu m}$ from the porous wall within the Upstream subdomain (see Fig. 8 below).

In both subdomains, the flow solution is initialized using constant states chosen so that at $t=0$, a shock occurs at $x=\Gamma$ with a specified initial speed $V_{s}$. Specifically:

• (a)

In the Upstream subdomain, the flow solution is initialized as follows

[TABLE]

where the notation $\spadesuit(V_{s},\textrm{RH})$ designates that $\spadesuit$ is determined from $V_{s}$ and the Ranking-Hugoniot jump conditions.

• (b)

In the Downstream subdomain, the flow is assumed to be initially quiescent and in a constant state determined by $V_{s}$ and the Rankine-Hugoniot jump conditions (see Tab. 1) – that is

[TABLE]

• (c)

The Mars atmosphere consisting mainly of carbon dioxide CO2, the gas properties are set to

[TABLE]

• 4.

At each subdomain boundary facing $\Gamma$ – that is, at each of the inlet and outlet boundaries – the boundary conditions are set to the uniform flow defined by the constant state used to specify the initial conditions in that subdomain.

• 5.

Two equivalent pore geometries defined by the aspect ratios $2/7\times 2/7$ and $8\%\times 100\%$ ($\alpha=8\%$) are considered.

• 6.

For each considered equivalent pore geometry, three unsteady numerical simulations are performed in a time-interval sufficiently large to capture the reflection of the shock wave on the porous wall:

• (a)

One three-dimensional pore-resolved simulation performed using DNS flow solver HAMeRS as described in Appendix B to generate reference data.

• (b)

Two counterpart one-dimensional simulations performed using the flow solver AERO-F [10, 8] equipped with the EBM FIVER [38, 19, 22, 12] and the homogenization approach presented in this paper (see Appendix A): one simulation performed using 400 grid points; and one performed using 2,000 grid points.

Each aforementioned supersonic instance of the porous flow model problem graphically depicted in Fig. 1 is defined by a specified shock speed value $V_{s}$. Three different values of $V_{s}$ are considered: $V_{s}=279.4\ \mathrm{m\ s^{-1}}$, $V_{s}=335.9\ \mathrm{m\ s^{-1}}$, and $V_{s}=367.1\ \mathrm{m\ s^{-1}}$. In all cases, the effect of bulk viscosity is neglected.

4.3.1 Pore-resolved numerical simulations

For each considered pore geometry and each considered shock speed value $V_{s}$, the results of the pore-resolved, unsteady simulation reveal that when the shock bounces off the porous wall, a high density/high pressure flow zone develops upstream of the porous wall. Consequently, the flow through the pore behaves like a flow through a converging-diverging nozzle. Downstream of the porous wall, the flow immediately accelerates and goes through an expansion. Fig. 9 displays the Mach number profiles computed at $t=23.46\ \mathrm{\mu s}$. It shows that as in the case of low-speed flow (see Fig. 6), the pore geometry with the smaller aspect ratio leads to jet streams. However, the reader can also observe that in the high-speed case, the flows are choked.

Let

[TABLE]

where $|\Omega^{\textrm{pore}}|$ denotes the area of the pore and $d=\sqrt{4|\Omega^{\textrm{pore}}|/\pi}$ is the pore diameter, denote the Reynolds number based on the jet flow velocity at the pore exit. Tab. 2 reports the values of this number computed at $t=23.46\ \mathrm{\mu s}$ by post-processing the results of the pore-resolved simulations. The reported low values, which are partially due to the low density and partially to the small pore size, highlight the importance of accounting for viscous effects when computing flows past porous walls

Fig. 10 reports for each considered pore geometry the variation of the averaged mass flow rate with the pressure difference between both sides of the porous wall computed by post-processing at each time-step of the computation the results of the corresponding pore-resolved, unsteady numerical simulation. Each variation is found to be linear and therefore describable by

[TABLE]

In other words, it satisfies Darcy’s law, where the constant $K$ depends on the pore geometry. In this case, a simple least square fitting yields $K=1.61\times 10^{-4}\ \mathrm{s\ m^{-1}}$ for the pore geometry with $2/7\times 2/7$ aspect ratio and $K=8.54\times 10^{-5}\ \mathrm{s\ m^{-1}}$ for that with $8\%\times 100\%$ aspect ratio.

4.3.2 Homogenization-based numerical simulations

Recall that all numerical simulations discussed herein are performed at the coupon level, due to the periodicity of the weave pattern (see Section 4.1). For this reason and due to the flat geometry of the considered porous wall (see Fig. 8), all counterparts of the pore-resolved, unsteady numerical simulations presented in Section 4.3.1 based on the homogenization approach described in this paper are performed in one dimension, in the computational fluid domain $x\in[-0.5,0.5]\ \mathrm{cm}$. Two discretizations of this domain in dual computational cells are considered for this purpose: one with $N=400$ grid points; and one with $N=2,000$ grid points. In each case, the porous wall is “embedded” in the discretization at $x=0$. The homogenization-based numerical simulations are performed using the vertex-based, finite volume flow solver AERO-F [10, 8] equipped with the second-order space-accurate EBM FIVER [22] (see also Appendix A) and the second-order explicit Runge-Kutta time-integrator constrained by the constant CFL number of 0.5.

Because for the equivalent pore with the aspect ratio of $2/7\times 2/7$ the jet stream was reported in Section 4.3.1 to expand immediately downstream of the porous wall (see Fig. 9), the thickness correction factor (see Section 3.3) is set for this case to $\eta_{f}=1$. For the equivalent pore with the aspect ratio of $8\%\times 100\%$, the jet stream can be seen in Fig. 9 to extend further in the computational fluid domain: thus, the thickness correction factor is set in this case to $\eta_{f}=4$.

Figs. 11, 12 and 13 report the averaged density, streamwise velocity, and pressure profiles computed at $t=23.46\ \mathrm{\mu s}$ using the proposed homogenization approach. It shows that: the one-dimensional discretization of $[-0.5,0.5]\ \mathrm{cm}$ in $N=400$ grid points is well resolved; and the results obtained using the proposed homogenization approach agree well with their counterparts obtained using the pore-resolved numerical simulations.

Next, attention is focused on the drag force $F_{D}$ generated by the porous wall. For a pore-resolved numerical simulation, $F_{D}$ can be defined as

[TABLE]

where $\Omega_{\textrm{vol}}$ denotes the volume representation of the porous wall, $\Omega_{\textrm{vol}}^{\textrm{pore}}$ denotes the part of $\Omega_{\textrm{vol}}$ occupied by the pores, $\partial{\Omega_{\textrm{vol}}}$ denotes the wet boundary surface of $\Omega_{\textrm{vol}}$, $\partial\Omega_{\textrm{vol}}^{\textrm{pore}}$ denotes the boundary surface of $\Omega_{\textrm{vol}}^{\textrm{pore}}$, and $n$ denotes the unit, outward normal to the surface over which the spatial integration is performed. Hence, $F_{D}$ consists in this case of a force due to the pressure jump across the porous wall (term ¡1¿ in (37)), a force induced by the jump of the normal component of the shear stress traction across the porous wall (term ¡2¿ in (37)), and the shear force on the pore boundaries (term ¡3¿ in (37)).

For a simulation based on the homogenization approach proposed in this paper, $F_{D}$ can be defined consistently with (37) as follows

[TABLE]

which in the case of the one-dimensional homogenized computations performed herein simplifies to

[TABLE]

Each of the integrals in (37) and (39) can be evaluated numerically using standard approximations. In the context of an EBM for CFD, the load evaluation algorithm presented in [19] can be used for this purpose.

Fig. 14 shows that for $N=400$, the homogenization-based numerical simulations accurately reproduce in all cases the variation with the shock speed of the drag force per unit area obtained at $t=23.46\ \mathrm{\mu s}$ using the counterpart pore-resolved numerical simulation.

5 Conclusions

The homogenization approach proposed in this paper for the treatment of porous wall boundary conditions distinguishes itself from alternatives in many ways. To begin, it is valid for both incompressible and compressible viscous flows. It does not require prescribing neither a mass flow rate nor a heuristic discharge coefficient. Instead, it models the inviscid flux through a porous wall surrounded by the flow as a weighted average of the inviscid flux at an impermeable surface and that through pores, and introduces a body force term in the governing equations to account for friction loss along the pore boundaries. The latter term is parameterized however by a thickness correction factor. The proposed approach, which can be implemented in any semi-discretization or discretization method, is successfully verified for a series of one-dimensional supersonic flow problems through a rigid porous wall, using three-dimensional pore-resolved numerical simulations. Its successful application to three-dimensional fluid-structure interaction problems associated with supersonic parachute inflation dynamics during Mars landing is reported in the companion papers  [14, 2].

Appendix A EBM FIVER and arbitrary embedding CFD meshes

Brief descriptions of the simplest version of the FInite Volume method with Exact two-material Riemann problems (FIVER) [38] and the implementation in this Embedded Boundary Method (EBM) for CFD of the proposed homogenization approach for the treatment of porous wall boundary conditions are provided here, in the context of arbitrary embedding CFD meshes. For more advanced versions of FIVER, the reader is referred to, for example, [22, 12, 4].

FIVER semi-discretizes the inviscid fluxes of Eq. 1 by a vertex-based Finite Volume (FV) method and the viscous fluxes of this equation by a FE method. Both methods are known for their ability to operate on arbitrary meshes. Using the standard characteristic function associated with a control volume (dual cell) $\mathcal{C}_{i}$ defined at a grid point $i$ of the structured or arbitrarily unstructured CFD mesh, the standard piecewise-linear tent test function $\psi_{i}$ associated with the grid point $i$, the standard bijection between the space of piecewise linear FE functions defined over the primal elements of the CFD mesh and the space of constant characteristic shape functions defined over the dual cells of this mesh [7], FIVER adopts the standard variational approach and integration by parts to semi-discretize Eq. 1 as follows

[TABLE]

In Eq. 40 above, $|\mathcal{C}_{i}|$ denotes the volume of $\mathcal{C}_{i}$, $\partial\mathcal{C}_{i}$ denotes its boundary, $\partial\mathcal{C}_{ij}$ denotes the cell facet shared by $\mathcal{C}_{i}$ and $\mathcal{C}_{j}$, $n_{ij}$ denotes the unit outward normal to $\partial\mathcal{C}_{ij}$, $\bm{W}_{h}$ denotes the spatial approximation of the fluid state vector $\bm{W}$, $\bm{W}_{i}$ denotes the average value of $\bm{W}_{h}$ in $\mathcal{C}_{i}$, $\mathcal{K}(i)$ denotes the set of grid points connected by an edge to the grid point $i$, $B$ denotes the computational fluid domain and $\partial B_{\infty}$ denotes its far-field boundary, $n_{\infty}$ denotes the unit outward normal to $\partial B_{\infty}$, $B^{e}$ denotes a primal element of the CFD mesh discretizing $B$, and $\psi^{e}_{i}$ denotes the restriction of $\psi_{i}$ to $B^{e}$.

Let $\Omega$ denote a surface embedded in the discretization of $B$ by a CFD mesh and let $(i,j)$ be an edge of this mesh that intersects $\Omega$ (see Fig. 15). In order to semi-discretize the inviscid flux through $\partial\mathcal{C}_{ij}$, the simplest version of FIVER assumes that the midpoint $M$ of the edge $(i,j)$ lies on $\Omega$ – and more specifically, $\partial\mathcal{C}_{ij}$ lies on $\Omega$; constructs a fluid-structure half-Riemann problem at the point $M$ in the direction of the normal $n_{ij}$ to $\partial\mathcal{C}_{ij}$; solves this half-Riemann problem to obtain the fluid state vector at the point $M$, $\bm{W}_{i}^{\star}$; and then approximates the inviscid flux through $\partial\mathcal{C}_{ij}$ as follows

[TABLE]

where $|\partial\mathcal{C}_{ij}|$ = $|\partial\mathcal{C}_{ji}|$ and $n_{ji}=-n_{ij}$. Note that by computing $\bm{W}_{i}^{\star}$ (or $\bm{W}_{j}^{\star}$) via the solution of a fluid-structure half-Riemann problem, FIVER is able to capture any shock or rarefaction wave near the embedded surface – that is, at the fluid-structure interface. Away from $\Omega$ – that is, when the edge $(i,j)$ does not intersect $\Omega$, FIVER approximates the inviscid flux through $\partial\mathcal{C}_{ij}$ using a standard FV method based on a standard approximate Riemann solver.

From (41) and Eq. 15, it follows that if $\Omega$ is a porous wall, the homogenized, numerical inviscid flux functions through $\Omega$ can be computed by FIVER as follows

[TABLE]

where averaging is performed as in (6). Again, the second-order extension of the above approximations using the classical Monotonic Upstream-centered Scheme for Conservation Laws (MUSCL) is straightforward.

As for the viscous approximation in Eq. 40, it can be re-written as

[TABLE]

due to the fact that the gradients of the piecewise-linear shape functions $\psi_{i}^{e}$ are constant in each primal element of the CFD mesh, $B^{e}$, whose volume is denoted here by $|B^{e}|$, and where $\bm{V}^{e}$ denotes an “element-value” of the vector of primitive fluid state variables. FIVER computes $\bm{V}^{e}$ and its gradient $\nabla\bm{V}^{e}$ as follows

[TABLE]

where $N^{e}$ denotes the number of grid points connected to the element $B^{e}$. If $B^{e}$ is intersected by an impermeable wall $\Omega$, FIVER reconstructs $\bm{V}^{e}$ by populating the velocity and temperature at the grid points connected to $B^{e}$ and occluded by $\Omega$ using constant or linear extrapolation procedures [19, 12]. This leads to the reconstructed vector $\bm{V}^{e,\textrm{be}}$. Hence, in general, FIVER computes $\bm{V}^{e}$ and $\nabla\bm{V}^{e}$ as follows

[TABLE]

It follows that if $\Omega$ is a porous wall, the homogenized approximation of the viscous fluxes at a grid point $i$ inside the porous wall can be computed as follows

[TABLE]

The body force term (see Eq. 30), which acts in the direction normal to the porous wall, is added at each grid point of the embedding CFD mesh as follows

[TABLE]

where $\bm{X}$ denotes position vector of a point in the computational fluid domain $B$, $n_{\textrm{w}}$ is the normal to the porous wall at the closest point $\bm{X}_{0}\in\Omega$ to $\bm{X}$,

[TABLE]

and $D(\bm{X})$ is the following Gaussian approximation of the Dirac function based on the distance $\|\bm{X}-\bm{X}_{0}\|_{2}$ to the porous wall

[TABLE]

where, as before, $\sigma=\eta_{f}/(2\pi)$.

Finally, it is noted that for the sake of computational efficiency, it suffices in practice to add the above source term only at the grid points located within one layer of elements of the embedding CFD mesh away from the porous wall.

Appendix B Pore-resolved numerical simulations

All pore-resolved numerical simulations discussed in Section 4 are performed using the flow solver HAMeRS [43, 42, 41] on uniform meshes. For each different simulation and pore geometry, the size of the computational domain graphically depicted in Fig. 1 – or Fig. 8, as appropriate – and its discretization are reported in Tab. 3. Specifically, each tabulated discretization is uniform in each direction and verified to be converged by comparing the results it delivers to their counterparts obtained on a discretization that is twice finer in each direction (for example, see Fig. 16).

Acknowledgments

Daniel Z. Huang and Charbel Farhat acknowledge partial support by the Jet Propulsion Laboratory (JPL) under Contract JPL-RSA No. 1590208, and partial support by the National Aeronautics and Space Administration (NASA) under Early Stage Innovations (ESI) Grant NASA-NNX17AD02G.