Modified symmetry technique for mitigation of flow leak near corners for compressible inviscid fluid flow

TL;DR

This paper introduces a modified symmetry technique to address flow leak issues near corners in compressible inviscid fluid flow simulations, improving accuracy in pressure prediction and shock structure modeling.

Contribution

The paper proposes a novel modification to the standard symmetry boundary condition technique to mitigate flow leaks near corners in Euler equation simulations.

Findings

Reduced flow leak near corners in numerical solutions.

Improved shock structure prediction on coarse meshes.

Enhanced accuracy in pressure estimation.

Abstract

Using the standard symmetry technique for applying boundary conditions for free slip and flat walls with corners will lead to flow leak through the wall near corners (violation of no penetration condition) and a corresponding error in prediction of pressure. Also, prescribing a state at the corner as a boundary condition is not possible. In this paper, a method for tackling the `corner point state' problem is given and modifications to the standard symmetry technique are proposed to mitigate flow leak near the corner. Using this modified symmetry technique, numerical solutions to the Euler equations for flows over forward facing and backward facing step are computed employing the Shu-Osher conservative finite difference scheme with WENO-NP3 reconstruction (with a formal order of accuracy in space of 3), Lax-Freidrichs flux splitting, and TVD-RK3 time discretisation. It is shown that using this modified symmetry technique leads to mitigation of flow leak near the corner and a better prediction of shock structure even on coarse meshes.

Tables

Table 1. Table 1 : Standard symmetry technique equations and corresponding boundary condition approximated at the wall boundary (see figure  3 for grid point labels). Subscripts are used to indicate grid point.

S.No.	SST Equations	Boundary condition that is approximated at point W
1	$p_{\tilde{K}}=p_K$	$\partial p/\partial x=0$
2	${\rho }_{\tilde{K}}={\rho }_K$	$\partial \rho /\partial x=0$
3	${\left(\rho u\right)}_{\tilde{K}}=-{\left(\rho u\right)}_K$	$\rho u=0$
4	${\left(\rho v\right)}_{\tilde{K}}={\left(\rho v\right)}_K$	$\partial \left(\rho v\right)/\partial x=0$

Table 2. Table 2 : L 1 subscript 𝐿 1 L_{1} and L ∞ subscript 𝐿 L_{\infty} errors of total energy density ( ρ ​ e t 𝜌 subscript 𝑒 𝑡 \rho e_{t} ) for different GPS for the problem of isentropic vortex moving in a uniform flow and observed order of accuracy.

GPS	$L_1\text{ error}$ $\times {10}^8$	$L_1\text{ order}$	$L_{\infty }\text{ error}$ $\times {10}^6$	$L_{\infty }\text{ order}$
1/25	862.397	-	1117.663	-
1/50	105.379	3.03	156.781	2.83
1/75	31.007	3.02	78.321	1.71
1/100	13.041	3.01	24.300	4.07
1/150	3.864	3.00	9.611	2.29
1/200	1.630	3.00	4.601	2.56
1/225	1.144	3.00	2.965	3.73

Table 3. Table 3 : Leakage in a Mach 4.0 flow over forward facing step: m ˙ l subscript ˙ 𝑚 𝑙 \dot{m}_{l} and e ˙ l subscript ˙ 𝑒 𝑙 \dot{e}_{l} as a percentage of m ˙ i subscript ˙ 𝑚 𝑖 \dot{m}_{i} and e ˙ i subscript ˙ 𝑒 𝑖 \dot{e}_{i} (Positive value indicates that mass or energy is flowing in).

GPS	Mass leak rate percentage (${\dot{m}}_l/{\dot{m}}_i\times 100$)			Energy leak rate percentage$({\dot{e}}_l/{\dot{e}}_i\times 100$)
	Below Corner ($a$)	Downstream of corner ($b$)	Total $\left(\left|a\right|+\left|b\right|\right)$	Below Corner ($c$)	Downstream of corner ($d$)	Total $\left(\left|c\right|+\left|d\right|\right)$
$1/50$	-0.09	0.36	0.45	-0.07	0.29	0.36
$1/100$	-0.04	0.20	0.24	-0.03	0.16	0.19
$1/200$	-0.01	0.09	0.10	-0.01	0.08	0.09

Table 4. Table 4 : Shock standoff distances for Mach 4.0 flow over a forward facing step for different WBCTs. SSTNGW predicts the grid independent shock standoff distance even for coarse meshes. ‘-’ indicates that computations for using a mesh with this GPS were not done (as grid independence was already achieved).

WBCT	Grid point spacing
	$1/50$	$1/100$	$1/200$	$1/400$	$1/800$	$1/1600$
SSTNGW	0.246	0.246	0.246	0.246	0.246	-
SSTGW	0.244	0.245	0.246	0.246	0.246	-
SSTNPE	0.246	0.247	0.246	0.246	0.246	-
SSTCVD	0.256	0.251	0.248	0.247	0.247	-
MST	0.281	0.264	0.255	0.251	0.249	0.248

Table 5. Table 5 : Mach 1.5 flow over backward facing step, mass and energy leak per unit time as a percentage of inflow mass and energy per unit time, respectively (Positive value indicates that mass or energy is flowing in).

GPS	Mass leak rate percentage (${\dot{m}}_l/{\dot{m}}_i\times 100$)			Energy leak rate percentage (${\dot{e}}_l/{\dot{e}}_i\times 100$)
	Below corner ($a$)	Upstream of corner ($b$) $\times {10}^3$	Total ($\left|a\right|+\left|b\right|$)	Below corner ($c$)	Upstream of corner ($d$) $\times {10}^3$	Total ($\left|c\right|+\left|d\right|$)
$1/50$	0.27	-3.3	0.28	0.37	-4.9	0.38
$1/100$	0.14	-1.8	0.14	0.19	-2.6	0.19
$1/200$	0.07	-0.8	0.07	0.09	-1.3	0.09

Equations

$$\frac{\partial }{\partial t}Q(x,t)+\frac{\partial }{\partial x}E(Q(x,t))=0$$

$$\frac{1}{\Delta x}\underset{-\frac{\Delta x}{2}}{\overset{\frac{\Delta x}{2}}{\int }}h(x+y,t)dy=E(x,t)$$

$$\frac{\partial E}{\partial x}\bigg{|}_{x=x_{o}}=\frac{h(x_{o}+\frac{\Delta x}{2},t)-h(x_{o}-\frac{\Delta x}{2},t)}{\Delta x}$$

$$\frac{\partial Q}{\partial t}\bigg{|}_{x=x_{0},t=t_{o}}+\frac{h(x_{o}+\frac{\Delta x}{2},t_{o})-h(x_{o}-\frac{\Delta x}{2},t_{o})}{\Delta x}=0$$

$$E^{\pm }=\frac{1}{2}\left(E(Q)\pm \alpha Q\right),$$

$$\alpha =\underset{Q}{\max }(|\vec{V}|+a),$$

$$\frac{\partial Q}{\partial t}\bigg{|}_{x=x_{0},t=t_{o}}+\frac{h^{+}(x_{o}+\frac{\Delta x}{2},t_{o})-h^{+}(x_{o}-\frac{\Delta x}{2},t_{o})}{\Delta x}+\frac{h^{-}(x_{o}+\frac{\Delta x}{2},t_{o})-h^{-}(x_{o}-\frac{\Delta x}{2},t_{o})}{\Delta x}=0,$$

$$\frac{1}{\Delta x}\underset{-\frac{\Delta x}{2}}{\overset{\frac{\Delta x}{2}}{\int }}{h}^{\pm }(x+y,t)dy=E^{\pm }(x,t).$$

$${\hat{h}}_{i+\frac{1}{2}}^{\pm }={\omega }_1^{\pm }{H}_1^{\pm }+{\omega }_2^{\pm }{H}_2^{\pm }$$

$$H_1^{+}=\frac{3E_i^{+}-E_{i-1}^{+}}{2},$$

$$H_1^{-}=\frac{E_i^{-}+E_{i+1}^{-}}{2},$$

$${\omega }_j^{\pm }=\frac{{\tilde{\omega }}_j^{\pm }}{{\bar{\omega }}^{\pm }},$$

$$\displaystyle\tilde{\omega}^{\pm}_{j}=\gamma^{\pm}_{j}\Bigg{(}1+\frac{\tau^{\pm}_{NP}}{(\beta^{\pm}_{j}+\epsilon)}\Bigg{)},~{}$$

$${\gamma }_1^{+}=\frac{1}{3},$$

$${\gamma }_1^{-}=\frac{2}{3},$$

$${\beta }_1^{+}=(E_{i-1}^{+}-E_i^{+}{)}^2,$$

$${\beta }_1^{-}=(E_{i+1}^{-}-E_i^{-}{)}^2,$$

$$\displaystyle\tau^{\pm}_{NP}=\Bigg{|}\beta^{\pm}_{3}-$$

$${\beta }_3^{+}=\frac{1}{4}(E_{i-1}^{+}-E_{i+1}^{+}{)}^2+$$

$${\beta }_3^{-}=\frac{1}{4}(E_i^{-}-E_{i+2}^{-}{)}^2+$$

$$E_k^{\pm }=E^{\pm }(x_k,t_n)$$

$$\frac{\partial }{\partial t}Q(x,y,t)+\frac{\partial }{\partial x}E(Q(x,y,t))+\frac{\partial }{\partial y}F(Q(x,y,t))=0,$$

$$\frac{d}{dt}u=L(u).$$

$$u^{n+1}=u^n+\Delta tL(u^n).$$

$$u^{(1)}$$

$$u^{(2)}$$

$$u^{n+1}$$

$$\frac{\partial Q}{\partial t}+\frac{\partial E}{\partial x}+\frac{\partial F}{\partial y}=0$$

$$Q=\left[\begin{array}{c}\rho \\ \rho u\\ \rho v\\ \rho e_t\end{array}\right],E=\left[\begin{array}{c}\rho u\\ \rho u^2+p\\ \rho uv\\ (\rho e_t+p)u\end{array}\right],F=\left[\begin{array}{c}\rho v\\ \rho vu\\ \rho v^2+p\\ (\rho e_t+p)v\end{array}\right],e_t=\frac{p}{\rho (\gamma -1)}+\frac{1}{2}\left(u^2+v^2\right)$$

$$u(x,y,t)$$

$$v(x,y,t)$$

$$\rho (x,y,t)$$

$$\frac{\partial }{\partial t}(\rho v)+\frac{\partial }{\partial y}(\rho v^2+p)+\frac{\partial }{\partial x}(\rho uv)=0\text{ (for CD)},$$

$$\frac{\partial }{\partial t}(\rho u)+\frac{\partial }{\partial x}(\rho u^2+p)+\frac{\partial }{\partial y}(\rho uv)=0\text{ (for BC)},$$

$$\frac{\partial }{\partial t}(\rho V_n)+\frac{\partial }{\partial n}(\rho V_n^2+p)+\frac{\partial }{\partial \tau }(\rho V_n{V}_{\tau })=0,$$

$$\frac{\partial }{\partial t}(\rho V_n)+\frac{\partial }{\partial n}(\rho V_n^2+p)+\frac{\partial }{\partial \tau }(\rho V_n{V}_{\tau })=0.$$

$$\frac{\partial }{\partial n}(\rho V_n^2)+\frac{\partial p}{\partial n}+\frac{\partial }{\partial \tau }(\rho V_n{V}_{\tau })=0,\text{ near corners}$$

$$\frac{\partial p}{\partial n}=\frac{\partial \rho }{\partial n}T+\rho \frac{\partial T}{\partial n}$$

Full text

1 Introduction

Problems with corners arise while computing numerical solutions to several differential equations like the heat equation, the incompressible Navier-Stokes equations and the compressible Euler equations. Applying boundary conditions at and near a corner is often an issue and it was addressed in many papers, some of which are:

• •

J. Crank et al [6] listed various strategies employed in tackling corners for linear elliptic and parabolic equations,

• •

H. Holstein. Jr [12], G. A. Ache [1], addressed the re-entrant corner issue for incompressible viscous fluid flow equations,

• •

Woodward et al [22] pointed out that the expansion corner of a free slip flat wall in the Mach 3.0 flow over a forward facing step is a singular point and they proposed a way for tackling problems that arise due to it,

• •

A. Verhoff [20] used analytic methods to address the problem of corners for compressible inviscid fluid flows, which were modelled using an approximation of the Euler equations.

In this paper, we revisit solutions to the compressible Euler equations. We describe problems with applying boundary conditions associated with free slip flat walls having corners that are not known to be stagnation points. We propose ways to address these problems.

Free slip condition corresponds to ‘no resistance to the flow at the boundary’. Solid and rigid wall boundary conditions are equivalent to a no penetration condition at the boundary. To enforce these conditions in computations, the wall is often assumed to be a surface of symmetry. Henceforth, the technique of enforcing the wall boundary to be a surface of symmetry will be referred to as the ‘standard symmetry technique’ or ‘SST’.

One issue with using SST is that it is accurate only for flat walls. Using SST leads to a zero pressure gradient normal to the wall ($\partial p/\partial n$), which is not always the case on a curved wall as pointed out by Moretti [17]. For example, in the two-dimensional flow over a wall with a radius of curvature $`R$’, $\partial p/\partial n=-(\rho V_{\tau}^{2})/R$ (where $V_{\tau}$ is the component of the velocity $\vec{V}$ that is tangential to the wall) is non-zero and therefore using SST will lead to inaccuracy in the numerical solution. The problem of applying an appropriate pressure gradient on the wall was addressed using different techniques, one amongst which is the ‘curvature corrected symmetry technique’, introduced by Dadone [7, 8].

Unfortunately, these techniques do not address the case of zero radius of curvature, namely a wall with a corner. A problem with the corner was addressed by Verhoff [20, 21], using analytic methods for solving approximations to steady, two-dimensional, compressible Euler equations, written in streamline coordinates. Approximate solutions for various problems including that of subsonic flow over a ramp [20], which has both compression and expansion corners, were reported. The first approximate solution obtained for the flow over a ramp predicted infinite momentum density at expansion corner and zero momentum density at compression corner. To tackle the singularity and obtain density, pressure, and magnitude of momentum, coordinate straining was used, which the author pointed out was ‘arbitrary to some extent’. Applying boundary conditions for such walls with corners is problematic for two reasons.

First, the normal, $\hat{n}$, and the tangent, $\hat{\tau}$, (see figure 1) at the corner are not defined. Consequently, the normal and tangential components of velocity, ($V_{n}$, $V_{\tau}$) and ($\partial p/\partial n$) can not be determined. Setting $\vec{V}=\vec{0}$ is not a solution, as the corner is not known to be a stagnation point. Simultaneous application of free slip and no penetration conditions does not seem possible at the corner. Further, in the pressure gradient equation ($\partial p/\partial n=-(\rho V_{\tau}^{2})/R$) the radius of curvature, $R$, is zero. Hence using these symmetry techniques, it is not possible to determine the state at the corner.

Second, using SST near a corner is inaccurate simply because, the wall, as will be shown in this paper, is not a surface of symmetry in the neighbourhood of the corner. A free slip wall boundary condition should preserve the condition of no penetration on the wall. In particular, the pressure on the wall should be such that it does not allow for flow penetration. It will be shown later that using SST will lead to flow leak through the wall near the corner. This implies that the pressure on the wall, in the neighborhood of the corner, obtained using SST, is erroneous.

In conclusion, there are two problems for free slip walls with corners, namely, tackling the ‘corner point state’ problem and ensuring no penetration on the wall in the neighborhood of the corner. These problems will be addressed in this paper.

Based on turning angle of the wall, we classify two-dimensional corners into compression and expansion corners as shown in figure 1. This study is restricted to Cartesian meshes and finite difference schemes. Consequently, turning angles are restricted to $\pi/2$ and $3\pi/2$ for compression and expansion corners, respectively. In case of a compression corner with turning angle of $\pi/2$, the velocity and the pressure gradient at the corner are zero. Therefore, there will be no error if SST is used near the corner. This leaves the problem of applying boundary conditions for expansion corners.

Two of the widely used test problems having walls with such expansion corners are, flow over a backward facing step and flow over a forward facing step. Numerical solution of flow over a backward facing step was published in a paper by Schmidt and Jameson [18]. Numerical solution of the Mach 3.0 flow over forward facing step (20% step height) was published in various papers, starting with Ashley F. Emery [9] and later in [14, 22, 13, 4]. Steady state solution of that problem was published by A. F. Emery and Bram Van Leer [9, 14], amongst others. In many of these papers, boundary conditions for walls with corners were applied using SST. For ease of demonstration of problems with using SST near corners, such as flow leak, we use the problem of Mach 4.0 flow over a forward facing step with 20% step height as this can be solved on a smaller domain, requiring less computational effort.

In this paper, we show that using SST for applying boundary conditions on walls near corners allows the violation of the no-penetration condition at and near the corner. We propose a method to tackle the ‘corner point state’ problem and modifications to the standard symmetry technique to ensure the no penetration condition is not violated near the corner. A solver employing the Shu-Osher conservative finite difference scheme with WENO-NP3 reconstruction, Lax-Freidrichs flux splitting and TVD-RK3 time discretisation is developed. This solver is used to compute supersonic flows over forward facing and backward facing step, using various boundary conditions near the corner and a comparison of the results obtained is presented here.

The rest of this paper is organized as follows. In section 2, the numerical method employed along with the standard symmetry technique is described. Using this scheme, the numerical solution of two test problems were computed for verification and validation, and the results obtained are presented. In section 3, the problem definition for flow over a forward facing step is given and mass leak is demonstrated. In section 4, the cause for the flow leak near expansion corners and ways to mitigate it are discussed. In section 5, modifications to the standard symmetry technique are proposed to reduce mass leak near corners, and in section 6, the modified symmetry technique is employed to obtain numerical solutions of Euler equations for flows over forward facing and backward facing step.

We begin with a brief description of the numerical method, the standard symmetry technique and their implementation.

2 Numerical method with the standard symmetry technique

In this section, the Shu-Osher conservative finite difference scheme with Weighted Essentially Non-oscillatory (WENO) reconstruction, Lax-Freidrichs flux splitting, and Total Variation Diminishing-Three stage Runge Kutta (TVD-RK3) time discretisation is described. This numerical method is used for solving all the problems presented in this paper. The standard symmetry technique used for simulating free slip flat walls is also described. We start with the description of the Shu-Osher conservative finite difference scheme.

2.1 Shu-Osher Conservative finite difference scheme

Consider a hyperbolic conservation law of the form

[TABLE]

Let the computational domain consist of grid points uniformly spaced in the physical domain, with grid point spacing equal to $\Delta x$. A function $h(x,t)$ is defined such that the sliding average of $h(x,t)$ over a length $\Delta x$ is equal to $E(x,t)$, that is,

[TABLE]

Taking a partial derivative of equation (2) with $x$, we get

[TABLE]

We refer to Barry Merriman [16] for detailed explanation and analysis of the Shu-Osher conservative finite difference scheme.

Using the method of lines and equations (1), and (3), a semi-discrete form of equation (1) is obtained at $x=x_{o},t=t_{o}$, which is

[TABLE]

The global Lax-Freidrichs flux splitting is described below.

2.2 Upwinding and Flux-Splitting

To account for propagation along the characteristic directions, upwind biasing of spatial derivatives is needed. This can be achieved by using flux splitting and appropriate biasing of the split fluxes. For flux splitting, we use the global Lax-Freidrichs flux splitting which is given below:

[TABLE]

where,

[TABLE]

where $a$ is the speed of sound and the maximum is taken over all the grid points in the computational domain. The semi-discrete form of the hyperbolic conservation law incorporating flux splitting becomes

[TABLE]

where

[TABLE]

The WENO-NP3 reconstruction procedure is used to obtain approximations to $h^{+}$ and $h^{-}$ using left and right biased stencils, respectively. It is described next.

2.3 WENO-NP3 reconstruction procedure

WENO-NP3 is one of the family of WENO reconstruction procedures. WENO reconstruction was introduced by Liu, Osher and Chan in 1994 [15]. Jiang et al gave a framework to build high order WENO schemes [13]. Changes to these schemes were proposed [2, 3, 10] to avoid loss of accuracy near critical points. One such scheme is the third order WENO-NP3 proposed by Wu et al [23], which maintains third order accuracy at critical points also. This will be briefly described below.

For WENO-NP3 reconstruction, a stencil of 3 points is used (see figure 2 for stencils and sub stencils). Equation (7) is used to advance from time $t_{n}$ to $t_{n+1}$. At grid point $x_{i}$, approximations $\hat{h}^{\pm}_{i+\frac{1}{2}}$ and $\hat{h}^{\pm}_{i-\frac{1}{2}}$ (subscript $n$, indicating time level, is dropped for brevity) to $h^{\pm}(x_{i+\frac{1}{2}},t_{n})$ and $h^{\pm}(x_{i-\frac{1}{2}},t_{n})$, respectively, are needed. These are given by the following equations:

[TABLE]

Formulae for $H^{\pm}_{1},H^{\pm}_{2},\omega^{\pm}_{1}\text{, and }\omega^{\pm}_{2}$ are given below,

[TABLE]

$H^{\pm}_{1}$, $H^{\pm}_{2}$ are lower order approximations to $h^{\pm}(x_{i+\frac{1}{2}},t_{n})$ and are calculated using the relevant sub-stencils shown in figure 2. The linear weights $\gamma^{\pm}_{1}$, $\gamma^{\pm}_{2}$ and the smoothness indicators $\beta^{\pm}_{1}$, $\beta^{\pm}_{2}$, and $\beta^{\pm}_{3}$ are used to calculate the nonlinear weights $\omega^{\pm}_{1}$, $\omega^{\pm}_{2}$. A convex combination of $H^{\pm}_{1}$, $H^{\pm}_{2}$, with the corresponding nonlinear weights is taken to obtain the final weighted essentially non-oscillatory reconstruction.

This procedure is for spatial discretisation of hyperbolic conservation law in one space dimension. For equations in two space dimensions such as,

[TABLE]

the same procedure can be used for discretising the $x$ and $y$ derivatives separately. The resulting semi-discrete form is integrated in time using TVD-RK3 method, which is described below.

2.4 TVD-RK3 time discretisation

Consider the equation

[TABLE]

The simple forward Euler time discretisation between two time levels $t_{n}$ and $t_{n+1}$ separated by $\Delta t$ is given by

[TABLE]

A three stage third order TVD (Total Variation Diminishing) or SSP (Strong Stability Preserving) [11] Runge-Kutta discretisation is given by

[TABLE]

The TVD-RK3 discretisation is used to advance in time from $t_{n}$ to $t_{n+1}$. Next, the standard symmetry technique used for applying wall boundary conditions is described.

2.5 Governing equations and the standard symmetry technique

We start with the two-dimensional Euler equations, which are

[TABLE]

where

[TABLE]

In finite difference schemes, wall boundary conditions can be applied by using ghost points (see Figure 3).

We classify meshes into two types, mesh with no grid points on the wall (MNGW - see figure 3(a)) and mesh with grid points on the wall (MGW - see figure 3(b)). For meshes with grid point on the wall, the state on the wall boundary is available either to directly apply the boundary condition, or, to verify if the applied boundary condition is producing the appropriate state on the wall. For a mesh without grid point on the wall, the boundary condition can be applied using the ghost points but the state at the boundary is not directly available and must be inferred.

The application of the standard symmetry technique at point W (see figure 3) involves setting the states at ghost grid points $\tilde{J}$ and $\tilde{K}$, which are located symmetrically with respect to the wall corresponding to grid points $J$ and $K$, respectively. The previous two sections have shown how the states at J, K and R are determined. The state at $\tilde{K}$ is found using the state at grid point $K$ and the equations in table 1 (subscripts of flow properties are used to indicate grid points).

The negation of normal component of momentum density in the ghost points ensures a zero normal component of velocity on the wall. The symmetry of tangential momentum will allow slip. The same process is used for other ghost points like $\tilde{J}$.

For code verification and validation, the test problems of isentropic vortex moving in uniform flow [19] and the problem of a shock reflection off a flat plate [24] are used.

2.6 Isentropic Vortex moving in uniform flow

The initial condition for this problem is an isentropic vortex perturbation added to a uniform flow in the positive $x$ direction. The solution at any time $t$ is given by:

[TABLE]

where $p(x,y,t)=(\rho(x,y,t))^{\gamma}$ and $r=\sqrt{(x-x_{0}-u_{0}t)^{2}+(y-y_{0})^{2}}$. The parameter values chosen are $x_{0}=8$, $y_{0}=0$, $\beta=2.0$, $u_{0}=1.0$ and $\gamma=1.4$. The computational domain is a square of dimensions $16\text{units}\times 16\text{units}$ with $0\leq x\leq 16$ and $-8\leq y\leq 8$. Periodic boundary conditions are applied along the $x$ and $y$ directions. Using the properties at $t=0$ as initial conditions, the numerical method described before is used to obtain a solution at $t=2.34375$ units. The time step ${\Delta t=0.09765625\Delta x}$. This problem was run for meshes with grid point spacings ($GPS=\Delta x=\Delta y$) of $1/25,1/50,1/75,1/100,1/150,1/175,1/200,\text{ and }1/225$. The $L_{1}$ and $L_{\infty}$ errors for meshes with different GPS and the observed order of accuracy are given in table 2.

2.7 Shock reflection off a flat plate

The problem of shock reflection off a flat plate [24] is used for validating SST and the numerical scheme. The governing differential equations are given by equation (28). The computational domain along with the boundary conditions are shown in figure 4. A shock, with shock wave angle of 29°, pre-shock Mach number of 2.9 reflects off a flat plate. The computational domain is a rectangle of size 3.5 units by 1.0 units in the $x$ and $y$ directions, respectively. The bottom boundary, ‘$y=0$’, is a free slip wall. The left boundary, ‘$x=0$’ is an inflow with $M=2.9$ flow in the positive $x$-direction, $p=1.0$, $\rho=1.4$. On and above the top boundary ‘$y=1.0$’, post oblique shock state is prescribed. The shock makes an angle of 151° with the $x$ axis, as shown in figure 4. The flow field is initialized with a Mach 2.9 flow with velocity vector pointed in the positive $x$-direction.

A CFL number of $0.125$ was chosen to calculate the global time step value. The ‘$\alpha$’ in the Lax-Freidrichs flux splitting is given by $\alpha=(2.9+1.0)=3.9$. Computations were done using MGW with grid point spacings or GPS ($=\Delta x=\Delta y$) of $1/50,1/100$, $1/200$, $1/400$, and $1/800$. Color plot of pressure obtained using MGW with GPS of $1/800$ is shown in figure 5.

The values of pressure, post the incident and reflected shock, of 2.14 and 4.10, obtained from the numerical solution are in good agreement with those of the analytic solution. The reflected shock in the numerical solution makes an angle of approximately 23° with the $x$-axis, which is also in good agreement with the analytic solution. Having tested the solver along with the standard symmetry technique, we turn to the problem of walls with corners.

3 Demonstration of flow leak near a corner

In this section, we demonstrate flow leak near the corner using the problem of Mach 4.0 flow over a forward facing step and we begin with the problem definition.

3.1 Mach 4.0 flow over a forward facing step: Problem domain and boundary conditions

Figure 6 shows the geometry of the flow field. Except for the inflow and outflow, all of the boundaries are free slip walls. In non-dimensional units $AG=1.0$, BC, the step height which should be $20\%$ of $AG$, is equal to $0.2$ and AB $=0.6$. The inflow is supersonic and the state of Mach 4.0 flow in the positive $x$ direction with $\rho=1.4$, $p=1.0$, is prescribed there. For ease of applying boundary conditions, CD $=0.6$ is chosen so that the outflow will be supersonic. One sided differences, biased in the negative $x$ direction are used to calculate $x$ derivatives near the outflow. The inflow conditions are prescribed at all grid points as the initial conditions. For those conditions, speed of sound $a=1$ and $\alpha\text{ (for Lax-Freidrichs flux splitting)}=5.0~{}(=a+M)$. As mentioned earlier, state at the expansion corner point ‘C’ (see figure 6) can not be determined and ways of tackling that problem are given next.

3.2 Tackling the ‘corner point state’ problem

As mentioned earlier determining the state at the corner is not possible.

One way to avoid this problem is to choose the mesh such that there is no grid point at the corner, but this results in the MNGW grid as shown in figure 7(a) and a corresponding approximate application of boundary condition.

To avoid this, we can have a mesh with grid points on the wall (MGW) as shown in figure 7(b) or 8, which will also have a grid point at the corner. Now, the problem of determining state at grid point ‘[math]’ needs to be addressed. Problems with corners whilst solving other partial differential equations were addressed in [12, 1, 20]. For the Euler equations on Cartesian mesh, we propose the following fix for the ‘corner point state problem’. In a Cartesian mesh, unlike compression corner points (like ‘$B$’ in figure 8), expansion corner points (like ‘[math]’ in figure 8) have all the necessary grid points required to discretise the Euler equations with an appropriate upwind biasing to the full order of the scheme. Therefore, we propose solving the governing differential equation at the corner, instead of applying boundary conditions. To repeat, no boundary condition is applied at ‘[math]’. This fix will be referred to as ‘corner fix’. The point ‘0’ becomes an interior point and the corner is, in a sense, ‘rounded’.

Using the corner fix implies that the no-penetration condition is violated in the portion of the boundary between grid points $-1$, [math] and also between [math], $1$ (see figure 8). This is because using the corner fix will lead to a finite non-zero velocity at the corner grid point [math]. Assigning any state with non-zero velocity at grid point [math] will lead to this violation. In addition to this, there will be flow leak in the region between grid points ‘$1$’ - $D$ and also between grid points ‘$-1$’ - $B$. This leak will be demonstrated using the solution for the Mach 4.0 flow over a forward facing step obtained using MGW and the ‘corner fix’, employing SST for applying wall boundary conditions.

3.3 Flow leak near corner

Numerical solution to the problem of Mach 4.0 flow over a step, described in section 3.1, is obtained using the numerical method described in sections 2.1 - 2.5.

The problem with the numerical solution obtained is that there is flow leak through the wall near the corner. Figure 9 has plot of magnitude of mass flux and energy flux through the wall boundary, below and downstream of the corner grid point. The grid points are labelled as shown in figure 8. As can be seen, the no-penetration condition is violated for several grid points below and downstream of the corner point and hence there is non-zero mass flux and energy flux through the wall boundary. This mass flux is a leak. It is integrated to find the total mass flow rate leaking, $\dot{m}_{l}$, near the corner. The total mass flow rate into the domain through the inlet boundary AG (figure 6), is denoted by $\dot{m}_{i}$. Similarly, $\dot{e}_{l}$ and $\dot{e}_{i}$ can be defined for the energy flow rate. The mass and energy leak data in table 3 is obtained by numerical integration (trapezoidal rule). The integration for leak below the corner is done from point ’B’ to grid point ’-1’ and that for downstream of corner is done from grid point ‘1’ to point ’D’ (see figure 8). This total mass (energy) leak per unit time is expressed as a percentage of inflow rate of mass (energy) through face AG (see figure 6), in table 3.

The flow leak reduces with reduction in grid point spacing which can be asserted using data in table 3. Positive values indicate that mass or energy is flowing into the domain. The flow leaks out of the domain below the corner and leaks into the domain downstream of the corner as is evident from the percentages in table 3.

Figure 10 shows plots of mass flux normal to the wall boundary near the corner for GPS of $1/50$, $1/100$ and $1/200$. Figure 11 has the corresponding plot of energy flux. From the plots, it is evident that the physical region over which flow leak happens decreases with decreasing grid point spacing.

The flow leak near the corner is due to the use of SST near the corner and this will be elaborated next.

4 Analysis of standard symmetry technique near corners

The analysis is performed using the grid points shown in figure 12. On the wall, the tangent and normal are defined except at grid point [math]. The limit of the equations for the normal component of momentum equation approaching the wall, for portions of the wall $CD$ and $BC$ are given by

[TABLE]

[TABLE]

respectively. For brevity, equations 33 and 34 will be represented using the following equation

[TABLE]

where $n$ (see figure 12) is measured in the direction normal to the wall, $\tau$ measured along the wall. Therefore while calculating approximation to $\tau$ derivatives at grid points on the line through $C$ and $D$, $(V_{n},V_{\tau})=(v,u)$. For $\tau$ derivatives at grid points on the line through $B$ and $C$, $(V_{n},V_{\tau})=(u,v)$.

Analytically, on portions of the wall BC and CD (see figure 12), it can be asserted that the term $\partial(\rho V_{n}V_{\tau})/\partial\tau$ is zero at every point except the corner C (grid point [math]). Also, as mentioned earlier, determining the state at [math] is problematic. In order to tackle this ‘corner point state problem’, we used the ‘corner fix’ (see section 3.2) of solving the discretised governing equations at [math], as is done at any interior grid point.

Let the velocity vector at grid point [math] be $\vec{V}_{0}$ and the normal at grid point $1$ be $\hat{n}_{1}$. Now, $\vec{V}_{0}.\hat{n}_{1}$ need not be zero, for a $\vec{V}_{0}$ obtained using the ‘corner fix’. In other words, velocity at grid point [math] need not be in the same direction as the tangent at $1$. Unfortunately, calculating an upwind biased approximation to $\partial(\rho V_{n}V_{\tau})/\partial\tau$ at grid point $1$, requires that the state at grid point [math] also be used. Approximation to $\partial(\rho V_{n}V_{\tau})/\partial\tau$, calculated using state at [math] will not be zero if $\vec{V}_{0}.\hat{n}_{1}$ is not equal to zero, which happens to be the case for the $\vec{V}_{0}$ obtained using the ‘corner fix’. Even if a different fix is used and a state at grid point [math] is assigned such that $\vec{V}_{0}.\hat{n}_{1}=0$, it will lead to a similar problem at grid point $-1$ unless $\vec{V}_{0}$ is set equal to $\vec{0}$. But grid point [math] is not a stagnation point. Therefore whatever non-zero velocity is assigned at grid point [math], it will lead to a non-zero $\partial(\rho V_{n}V_{\tau})/\partial\tau$ at either grid point $1$ or $-1$ or at both of them.

Using SST at grid point $1$ will lead to $\partial(\rho V_{n}^{2}+p)/\partial n=0$. $\partial(\rho V_{n}V_{\tau})/\partial\tau\neq 0$ and $\partial(\rho V_{n}^{2}+p)/\partial n=0$ will lead to a nonzero normal component of momentum at grid point $1$ through equation (35). Repeating this argument, this non zero normal component of velocity at grid point $1$ will result in a non zero normal component of velocity at grid point $2$. In subsequent time steps, this nonzero normal velocity will propagate to other grid points along the wall. The same will happen at grid points $-1$, $-2$, …, which are below grid point [math]. Different ways to address this flow leak for these problems are described next.

4.1 Algorithmic fixes for the flow leak problem

If no-penetration condition is satisfied at grid points $-1$ and $1$, the nonlinear WENO weights will make sure that $\partial(\rho V_{n}V_{\tau})/\partial\tau$ is essentially zero at grid points below $-1$ and grid points to the right of $1$. This, along with SST will prevent mass leak along the wall at all grid points except the corner. Two algorithmic fixes to achieve this are given next.

The first fix is to set the normal component of velocity to zero (enforcing no-penetration), at grid points $-1$ and $1$, after each time step or Runge-Kutta stage. This, as mentioned earlier, will prevent the flow leak from happening at grid points below $-1$ and grid points to the right of $1$. This technique will be referred to as ‘SSTNPE’ (Standard symmetry technique with no penetration enforced).

The second fix is called the corner velocity direction fix. Let the density, momentum density, and total energy density at the grid point [math], obtained using the ‘corner fix’ be $\rho_{0},\,\rho_{0}\vec{V}_{0},\,\rho_{0}e_{t_{0}}$, respectively. Let $|\vec{V}_{0}|=V_{0}$. Let $\hat{e}_{x}$ and $\hat{e}_{y}$ be the unit vectors along the positive $x$ and $y$ directions, respectively. There are two grid lines through grid point [math] as shown in figure 13, one parallel to $x$-axis (figure 13(b)) and one parallel to $y$-axis (figure 13(a)). For calculating $x$ and $y$ derivatives on these grid lines, except at grid point [math], the following state at grid point [math] is used: $\rho_{0},\,\rho_{0}V_{0}\hat{\tau_{0}},\,\rho_{0}e_{t_{0}}$. As shown in the figure 13, $\hat{\tau_{0}}=\hat{e}_{y}$ for grid points along the grid line parallel to $y$-axis and $\hat{\tau_{0}}=\hat{e}_{x}$ for grid points along the grid line parallel to $x$-axis. That is, while calculating $x$ derivatives (like $\partial(\rho u^{2}+p)/\partial x$ and $\partial(\rho uv)/\partial x$) at grid points like $H$ and $1$ (figure 12), the velocity at grid point [math] is taken to be $V_{0}\hat{e}_{x}$. Similarly, while calculating $y$ derivatives (like $\partial(\rho v^{2}+p)/\partial y$ and $\partial(\rho uv)/\partial y$) at grid points like $I$ and $-1$, the velocity at grid point [math] is taken to be $V_{0}\hat{e}_{y}$. The WENO weights and prescribed direction of velocity will ensure $\partial(\rho V_{n}V_{\tau})/\partial\tau$ will be zero at grid points $-1$ and $1$, which along with SST will ensure no-penetration condition is satisfied at grid points $-1$ and $1$. This fix will be referred to as ‘SSTCVD’ (Standard symmetry technique with corner velocity direction fix). This is similar to the suggestion of Verhoff [20] that “At the corner points the velocity (or momentum) vector rotates at constant magnitude through an angle $\alpha$ [$\alpha$ is the flow turning angle which is equal to $\pi/2$ for problems considered in this paper] due to an impulsive-type interaction”.

Using either ‘SSTNPE’ or ‘SSTCVD’, the no-penetration condition will be satisfied at all grid points on the wall except the corner. However, it will lead to contribution of the non-zero term ‘$\partial(\rho V_{n}V_{\tau})/\partial\tau$’, at grid points $1$ and $-1$, being ignored and a corresponding error in the state at grid points $1$ and $-1$. To avoid this, we retain this non-zero derivative and propose a new technique to apply wall boundary conditions at grid points $1$ and $-1$. In the next section, we propose modifications to the standard symmetry technique so that the effect of the non-zero term, $\partial(\rho V_{n}V_{\tau})/\partial\tau$, is also considered.

5 Modified symmetry technique for walls with expansion corners

Now, we propose modifications to the standard symmetry technique to incorporate the non-zero tangential derivative, $\partial(\rho V_{n}V_{\tau})/\partial\tau$, near corners and derive equations to be used at grid points near corners for applying boundary conditions.

5.1 Condition on normal derivative of pressure on the wall, near corners

The governing differential equations are solved at the corner grid point (point 0 in figure 12). The equations that will be derived next, are for applying boundary conditions at grid points adjacent to the corner, that is for grid points $-1$ and $1$. At grid points near the corner, the normal direction is defined. The limit of the normal momentum equation approaching the wall is given by equation (35) and is repeated below.

[TABLE]

Near the corner, on the wall, no penetration implies $\rho V_{n}=0$ and $\partial(\rho V_{n})/\partial t=0$. That leaves us with an equation

[TABLE]

To ensure that the normal component of velocity on the wall is zero, the states in the ghost points must be such that they satisfy the discretised version of equation (36). It is pointed out that while discretising equation (36), the derivatives should be calculated using the Lax-Freidrichs flux splitting.

For a third order scheme, there will be 2 layers of ghost points (as shown in figure 12) and values for density, pressure, and velocity are needed at these points in these 2 layers. For ensuring free slip we use the equations in the standard symmetry technique written in section 2.5 except for pressure and density. As mentioned earlier, discretised form of equation (36) is used to calculate pressure gradient. We also need the gradient of density (or temperature) to define the states at the ghost points. A zero normal temperature gradient is chosen for applying the boundary conditions. To calculate density gradient normal to the wall the equation of state $p=\rho RT$, with $R=1$ is used.

Taking a derivative of the equation of state $p=\rho T$ along $\hat{n}$, we get

[TABLE]

Setting gradient of temperature normal to wall in equation (37) to zero, we get the following:

[TABLE]

Using the discrete form of equations (36) and (38) at $-1$, we must determine pressure and density at grid points at $\tilde{J}$ and $\tilde{K}$. This is not possible as there are only 2 equations but 4 unknowns. To avoid this problem we can set pressure and density at one of the points $\tilde{J}$ and $\tilde{K}$ equal to that of $J$ and $K$, respectively, that is $p_{\tilde{J}}=p_{J},\rho_{\tilde{J}}=\rho_{J}$ or, $p_{\tilde{K}}=p_{K},\rho_{\tilde{K}}=\rho_{K}$. Choosing $p_{\tilde{J}}=p_{J},\rho_{\tilde{J}}=\rho_{J}$ is better because depending on the WENO weights, the normal pressure derivative on the wall may be independent of pressure and density at grid point $\tilde{J}$. Therefore, we choose the equation $p_{\tilde{J}}=p_{J},\rho_{\tilde{J}}=\rho_{J}$ to eliminate the pressure and density at $\tilde{J}$ which leaves us with two variables and two equations which can be solved. Using discrete forms of equations (36) and (38), we solve for pressure and density at the grid point $\tilde{K}$. A linear central difference with formal order of accuracy of 4, was used to discretise equation (38) and the Shu-Osher conservative finite difference method with WENO-NP3 reconstruction and Lax-Freidrichs flux splitting was used to discretise equation (36). For WENO-NP3 reconstruction, equation (36), is non-linear in $\rho_{\tilde{K}}$ and $p_{\tilde{K}}$. Bisection method was used to solve discretized forms of equations (36) and (38) to an accuracy of $10^{-8}$. This modified symmetry technique described above will be referred to as MST.

5.2 Algorithm for solving for pressure in MST using bisection method [5]

Choose an initial bounding interval for $p_{\tilde{K}}$: $[a,b]\text{, where }a<b$ 2. 2.

Calculate densities at ghost point $\tilde{K}$ ($\rho_{\tilde{K}}$) using equation (38) with $p_{\tilde{K}}=a$ and $p_{\tilde{K}}=b$ and let them be $\rho^{a}_{\tilde{K}}$, $\rho^{b}_{\tilde{K}}$, respectively. 3. 3.

Check if setting $p_{\tilde{K}}=a,\rho_{\tilde{K}}=\rho^{a}_{\tilde{K}}$ produces numerical approximation to LHS of equation (36) with opposite sign to that obtained by setting $p_{\tilde{K}}=b\text{ and }\rho_{\tilde{K}}=\rho^{b}_{\tilde{K}}$. Else, go back to step 1 and change the bounding interval. 4. 4.

Using the bounding interval of pressure ($[a,b]$) and using a 4th order accurate discretisation of equation (38), calculate three densities at $\tilde{K}$ corresponding to $p_{\tilde{K}}=a$, $p_{\tilde{K}}=(a+b)/2$ and $p_{\tilde{K}}=b$ and let those densities be $\rho^{a}_{\tilde{K}},\rho^{(a+b)/2}_{\tilde{K}},\rho^{b}_{\tilde{K}}$, respectively. 5. 5.

Calculate three numerical approximations to the LHS in equation (36) using flux splitting according to the procedure given in sections 2.1 - 2.3, by using ${(p_{\tilde{K}},\rho_{\tilde{K}})=(a,\rho^{a}_{\tilde{K}})}$, ${(p_{\tilde{K}},\rho_{\tilde{K}})=((a+b)/2,\rho^{(a+b)/2}_{\tilde{K}})}$, ${(p_{\tilde{K}},\rho_{\tilde{K}})=(b,\rho^{b}_{\tilde{K}})}$ and label them $R_{1},R_{2},R_{3}$, respectively. 6. 6.

Using $R_{1},R_{2}\text{, and }R_{3}$, choose a new, smaller bounding interval for $p_{\tilde{K}}$. If $(R_{1}R_{2})<0$ the new interval is $[a,(a+b)/2]$, otherwise it is $[(a+b)/2,b]$. 7. 7.

With the new interval obtained in the above step, repeat from step 4 till values of pressure and density are obtained such that equation (36) is satisfied with an error of $10^{-8}$.

The initial bounding interval for pressure for starting the bisection method is chosen as the interval $[-\eta p_{prev},\eta p_{prev}]$, where $\eta>0$ and $p_{prev}$ is the root of equation (36) in previous time step or RK Stage and $\eta$ is chosen so as to satisfy the condition in step 3.

5.3 A note on flow field initialisation and boundary conditions and CFL number

For ease of initialisation, the flow field may be initialised with uniform flow at all grid points including the grid points on the wall. This will lead to the normal component of velocity not being zero at some grid points on the wall, initially. In such a case, in addition to applying MST, the normal component of velocity must be set to zero after every RK stage for a few thousand time steps until the normal component of momentum on the wall settles to zero or a very low value. Also, for the first few thousand time steps, a CFL number of $1/128$ should be used and later (after 5000 or 10000 time steps) it can be increased to a higher value like $1/8$ or $1/4$. Failing to do this may lead to severe convergence problems.

The algorithmic fixes and modified symmetry technique described in previous sections are used to obtain numerical solution to flows over backward and forward facing step. A comparison of solutions obtained using these different techniques is presented in the next section.

6 Testing the new boundary technique (MST) and algorithmic fixes

The modified symmetry technique (MST) described in section 5 and the algorithmic fixes described in section 4.1 are tested by solving supersonic flows over forward facing and backward facing step. A comparison of numerical solutions obtained using these different techniques for meshes with different grid point spacings is presented. Labels for the five different wall boundary condition techniques (WBCTs) are given below:

1. (a)

Standard symmetry technique with no grid points on the wall and at corner - SSTNGW, 2. (b)

Standard symmetry technique with grid points on the wall and corner, and the corner fix being used - SSTGW, 3. (c)

Standard symmetry technique with grid points on the wall, with the corner fix being used, and no penetration enforced near the corner - SSTNPE (refer to section 4.1), 4. (d)

Standard symmetry technique with grid points on the wall, with corner fix being used, and modification of corner velocity direction - SSTCVD (refer to section 4.1), 5. (e)

Modified symmetry technique - MST (refer to section 5).

We start with the flow over forward facing step.

6.1 Mach 4.0 flow over a forward facing step

The problem domain and boundary conditions were described in section 3.1.

Figure 14 has colour plots of density obtained using meshes with GPS of $1/800$ and all the five WBCTs mentioned in the previous section. All the five techniques result in a Mach reflection on the top wall (GE, see figure 6). Similar solutions are obtained for GPS of $1/400$. Figure 15 has the colour plots of density obtained using meshes with GPS of $1/200$. Apart from MST and SSTCVD, the other techniques fail to produce a Mach reflection on the top wall, as they did for meshes with GPS of $1/400$ and $1/800$. Using MST resulted in a Mach reflection on the top wall for all tried GPS of $1/50,1/100,1/200,1/400,\text{ and }1/800$, whereas using SSTCVD did not produce a Mach reflection for GPS of $1/50$ and $1/100$, as shown in figure 16.

Table 4 has the shock standoff distances (SB in figure 6) for the 5 WBCTs, for different grid point spacings. The portion of shock near the wall AB is a normal shock with pre and post shock densities of $1.400$ units and $6.474$ units, respectively. The point S (shock location, in figure 6) is taken to be located on AB where the density is equal to half of the pre and post shock densities (which is equal to $3.937$). For SSTNGW, SSTGW, SSTNPE and SSTCVD, the shock standoff distance is grid independent for GPS of $1/400$ and $1/800$. The same is true for MST for GPS of $1/800$ and $1/1600$. With this, the grid independence is achieved for all the WBCTs used.

It is evident from table 4 that SSTNGW predicts the grid independent shock standoff distance even for coarser meshes (with GPS of $1/50$). However, as shown in figures 14 and 15, except for MST and SSTCVD, none of the other WBCTs are able to accurately capture the grid independent shock structure - that of a Mach reflection on GE (see figure 6). MST captures this shock structure even for GPS of $1/50$ and $1/100$, whereas SSTCVD does not, as shown in figure 16.

Next, the WBCTs are used for computing flow over backward facing step.

6.2 Supersonic flow over a backward facing step

Figure 17 shows the sketch of problem domain and boundary conditions. The lengths of different portions of the flow field are (refer to figure 17 for labels) - AB = $0.2$ units, AG = $0.8$ units, BC = $0.2$ units, ED = $1.0$ units. Except for the inflow and outflow, all of the boundaries are free slip walls. The inflow is supersonic. For ease of applying boundary conditions, the length of CD is chosen so that the outflow is also supersonic. In our computations, for different WBCTs, different values for CD in the range of $3.1$ units to $4.8$ units were chosen. At the inflow boundary, the state is prescribed. At the outflow boundary one sided differences, biased in the negative $x$ direction are used to calculate $x$ derivatives.

Conditions at inflow are: $\rho=1.4$, $(u,v)=(M,0)$, where $M$ is the inflow Mach number, and $p=1.0$. Computations for two inflow conditions with Mach numbers of 1.5 and 2.5 were done using meshes with GPS of $1/50,1/100,1/200,\text{ and }1/400$ for the five different WBCTs.

For both inflow Mach numbers of $1.5$ and $2.5$, the five WBCTs produce similar numerical solutions all having an expansion fan and reattachment shock, which reflects off the top wall. The solutions for inflow Mach number $1.5$ differ in the position of the reattachment shock for different WBCT. The same trend appears for the inflow Mach number of 2.5.

Next, the flow leak due to using SST for this problem is described and it is compared with that of the Mach 4.0 flow over a forward facing step.

6.2.1 Flow Leak

Table 5 has the mass and energy leak near the corner as a percentage of the inflow mass and energy, for the Mach 1.5 flow over a backward facing step for grid point spacings of $1/50$, $1/100$ and $1/200$. These are calculated similar to the data in table 3 using trapezoidal rule. The leak percentages are similar for the Mach 2.5 flow also. Figure 18 has plots of mass flux on the wall boundary near the corner for the Mach 1.5 and 2.5 flows for comparison.

The mass and energy leak percentages for flow over backward facing step are similar to that of the Mach 4.0 flow over forward facing step, as can be seen from the data in tables 5 and 3. A major portion of the leak happens below the corner for the flow over a backward facing step, whereas for flow over forward facing step it happens downstream of the corner. For Mach 4.0 flow over forward facing step, the leak below the corner is approximately one order of magnitude less than the leak downstream of the corner (see table 3), whereas for Mach 1.5 flow over backward facing step, the leak upstream of the corner is approximately two orders of magnitude less than that below the corner.

7 Conclusions

The problems of state at the expansion corner point and flow leak due to using standard symmetry technique near corners were addressed. A method to tackle the ‘corner point state’ problem was proposed. Using MGW, it was shown that using SST (SSTGW) will lead to leak near the expansion corner (see figures 10 and 11) and that refining the mesh will lead to reduction of flow leak near the corner as evidenced by the data in table 3.

To reduce the flow leak and to limit it to the corner point, three WBCTs - SSTNPE, SSTCVD and modifications to the standard symmetry technique (MST) - were proposed and implemented. The problem of leak at the corner still exists in the new WBCTs proposed. It is not clear how it can be eliminated because the normal and tangent at the corner are not defined and simultaneous application of free slip and no-penetration at the corner is not possible.

Results obtained using the five different WBCTs for flows over forward facing and backward facing step were presented and compared. Of the five WBCTs, for SSTNPE, SSTCVD and MST there is no mass leak at any grid point on the wall except the one at the expansion corner. Of SSTNPE, SSTCVD, and MST, only MST takes into account the term $\partial(\rho V_{n}V_{\tau})/\partial\tau$ (in equation (35)) for enforcing no-penetration condition, while SSTNPE and SSTCVD do not.

For the Mach 4.0 flow over forward facing step, SSTNGW predicts the grid independent shock standoff distance for the coarser meshes also. MST gives a better prediction of the shock structure (the type of shock reflection that happens at the wall GE, see figure 6). Using MST, a Mach reflection at the wall GE (see figure 6) and a $\lambda$ shock was obtained for all grid point spacings used (see figure 16). Whereas for the other corner techniques, only the finer meshes gave a solution with Mach reflection and $\lambda$ shock (see figure 14). The solutions obtained with coarser meshes have regular shock reflection (see figure 15, 16).

For the problem of flow over a backward facing step, the total mass and energy leak percentages due to using SST were similar to that for flow over forward facing step. For flow over backward facing step the major portion of the leak happens below the corner whereas this happens downstream of the corner for flow over forward facing step.