The issue on boundary conditions for hyperbolic problems and particularly for compressible fluid flows is a classic topic and so is the corresponding numerical treatment. There are a number of contributions via various approaches in literature, which roughly consist of three types of concerns: physical considerations, mathematical justifications (well-posedness arguments)  and numerical treatment. Physical considerations prescribe boundary data for a part of physical variables based on specific problems such as the solid-wall boundary condition; mathematical well-posedness justifies the validity of modelings subject to the prescribed boundary conditions; while numerical boundary conditions are prescribed for all physical variables so that discrete (approximate) equations can be implemented practically. These concerns, though with different objectives, have the common goal on correctly describing the underlying problems, for which the compatibility among the governing equations, prescribed boundary conditions and the initial data is a fundamental issue. As far as the numerical treatment is concerned, extrapolation technique is often employed, particularly for high order accurate numerical methods. For example, in [17, 18]
a lagrangian interpolation is performed to achieve a second order accurate approximation to the boundary data in space. However, it just gives the first order accurate approximation in time. Other works can be found, e.g. in[7, 15], in the finite volume framework, and even in complex geometries . In , characteristic method is used but restricted to first order accuracy for smooth flows. Although the resulting schemes may be well-implementable, the corresponding validation is not clear both from rigorous mathematical analysis and numerical performance. Improper extrapolation may lead to numerical instability such as spurious oscillations . From the viewpoint of numerical analysis, it is questionable whether the numerical boundary conditions are compatible with the discretized governing equations even though the underlying PDE models are well-posed. Hence it is worth addressing issue even though there are lots of studies available [7, 11, 12, 13, 14, 15, 16, 17, 18, 19, 28, 29, 30].
and focus on the left boundary . We assume that for example and inspect various situations upon the boundary requirement on . There are three typical cases:
, . For this case, the solution contains a rarefaction wave
, . For this case, we have a shock solution,
, . There exists no physically admissible solution for such a case and so the boundary condition is not well prescribed.
This example shows the subtlety of nonlinear problems as investigated in , unlike linear hyperbolic problems. In fact, for linear problems, physical boundary conditions depend on characteristic propagations. While for nonlinear compressible fluid flows, many physical boundary conditions are prescribed upon surroundings and cannot be even given a priori, such as the interaction of shock with solid boundaries  and the solid body floating in the air . Numerically, situations become more complicated. First, numerical boundary conditions should be given for all variables in order to be suitable for the computation so that proper extrapolations have to be used. For the strong interaction of waves with physical boundaries, the nonlinearity actually prevents the validity of extrapolations that may result in factitious phenomena. Second, high order approximations of boundary conditions are often made independently of the discretization of the governing equations, which may lead to incompatibility and loss of accuracy.
The GRP formulated here is different from the traditional GRP [3, 4] and more suitably called one-sided GRP. It is an initial-boundary value problem rather than a purely initial value problem. Correspondingly, a numerical method to solve this problem is called a one-sided GRP solver. Such a study has two-fold goals: (i) It is used to the compatibility of prescribed boundary condition with the governing equations; (ii) It proposes a family of high order numerical boundary conditions effectively compatible with the discretized governing equations. In fact, accelerated piston problems [8, 13] are the one-sided GRP formulated here and they have been refined to put into the simulation of fluid flows with moving boundaries [14, 13]. In , the similar idea was employed if the flow is smooth and consistent with the inverse Lax-Wendroff method . Note that since the two-stage fourth order framework [23, 24, 12] can be used to develop high order methods, there is no need to compute derivatives of order more than second. Hence this paper just focuses on second order GRP solvers, which provides a reliable tool no matter whether the solution is smooth or not.
We organize this paper as follows. In Section 2, we formulate the one-sided GRP problem. In Section 3 we discuss numerical boundary conditions for compressible Euler equations via the passage of one-sided GRP. We implement the resulting scheme in Section 4 and particularly display numerical results in Subsection 4.2 to demonstrate the performance.
2. High Order Numerical boundary conditions and one-sided GRP solver for hyperbolic balance laws
In this section we formulate one-sided Riemann problems for hyperbolic balance laws in a general framework and discuss related numerical solvers for the construction of numerical boundary conditions. This is different from the classical initial-boundary value problem for compressible fluid flows in  since initial and boundary values are generally not compatible in a continuous way. Such a setting is proposed for practical request. For example, the flow is discontinuous at the reflection point of a shock on a solid boundary. Correspondingly, the associated initial-boundary value problem is formulated below as the one-sided generalized Riemann problem (OS-GRP).
Consider hyperbolic balance laws
where is a source term representing external forces or geometrical effects, is the flux function. This system includes the compressible Euler equations we specified in the next section and many other models . It is assumed to be hyperbolic in the sense that the Jacobian of has
with a complete set of associated eigenvectors,
Each is genuinely nonlinear or linearly degenerate in the sense of Lax .
We focus on the left boundary . The right boundary is treated similarly. We emphasize that the free boundary problem can be studied too . On the boundary , the data is imposed as
where the operator projects the solution onto the boundary , is the number of negative eigenvalues. If (2.1) is a linear problem, i.e., , is a constant matrix, then the operator can be expressed in the matrix form
if has positive eigenvalues, where is a matrix. For nonlinear problems the eigenvalues depend on the solution and thus the integer may vary depending on the solution too. Hence the precise meaning of the operator is determined together with the solution of (2.1), as shown for the Burgers equation.
Denote the trace operator on the boundary ,
Then we propose the following assumption.
in some “appropriate” sense.
This assumption is very “rough” and understood in certain intuitive way. A primary judgement of the well-posedness boils down to the solvability of the following one-sided Riemann problem, while the dynamics is dependent on the one-sided generalized Riemann problem (GRP). Numerically each component of should be given a value on the boundary so that the corresponding numerical code can be implemented. Even with extrapolation, the approximation should be consistent with the solution of this one-sided GRP up to some desired accuracy order.
In this section, we denote by the boundary value for the solution , and the derivative of along the boundary .
2.1. Linear equations with constant coefficients
Let’s first get motivation from linear equations. Consider with the assumption as above
The characteristic decomposition tells that
where , is a left-eigenvector associated with the eigenvalue . The solution formula is
To obtain the solution values on the boundary , we have from (2.8)
Hence the boundary value of can be obtained by solving the following system
, are the row vectors of the matrix. Such a solution formula in turn helps to develop high order schemes. We can refer to [12, 29] and next sections for the practical implementation in gas dynamics, corresponding to the acoustic case of one-sided GRP problem.
The above discussion is of course made in the theoretical viewpoint. Numerically, we implement at each time level , as follows.
First order approximation. The initial data is assumed to be constant . Then we derive all components of by solving the following system
Obviously, this is exactly the same as the usual extension from the neighboring interior point using the characteristic method.
Second order approximation. As high order approximations are concerned, we not only need to know the value in the first order approximation, but we have to approximate the value as well. We denote by and subsequently . Then we have
Solving this system yields the value . This second order approximation shows clearly that the source term is input into the numerical boundary condition, unlike some direct extrapolation technique. Moreover, this characteristic method allows to deal with discontinuities at the origin.
Indeed, these two approximations correspond to the one-sided Riemann problem and one-sided generalized Riemann problem, respectively.
2.2. One-sided Riemann problem
The boundary data is prescribed as
In order to solve this problem, we can mimick the method for the standard Riemann problem in the state space [9, 20, 31]. At least for Euler equations, we will show how to solve it in the next section. The solvability of this one-sided Riemann problem depends on the compatibility of the prescribed boundary conditions with the initial data. Generally speaking, as shown for the Burgers equation, this problem may not have to be well-posed. Hence this one-sided Riemann problem plays a role in checking whether the boundary conditions are correctly prescribed.
Another role of the one-sided Riemann problem is to supplement all state variables for the practical calculation because the boundary conditions just prescribe partial components of them. For instance, consider the linear case, as indicated in (2.4), with characteristics leaving the boundary so that the rank of the boundary operator is . Then we use the characteristic decomposition to obtain other equations, as shown above.
Assume that we are able to solve this problem and obtain the solution with the trace on the boundary such that,
Then it is necessary to check whether there are exactly characteristics leaving from the boundary , similar to the linear case.
That is, the dimension of manifold Just like the case for the Burgers equation, this is not necessary true. Hence the solvability of one-sided Riemann problem is a necessary to judge the well-posedness of initial-boundary value problem for (2.4).
2.3. One-sided GRP
As (2.1) includes a source term or/and the initial condition is not uniform (typically consists of piecewise polynomials), one has to consider a one-sided generalized Riemann problem (GRP). From the numerical point of view, one needs to have high order accurate prescription of all components of on the boundary as well as the construction of spatial variation near the boundary when high order methods are sought.
For completeness, the one-sided generalized Riemann problem (GRP) is reformulated here as the initial and boundary value problem,
where is smooth, and is measurable. This is associated with the one-sided Riemann problem above. Similar to the interrelation between the standard generalized Riemann problem and the associated Riemann problem, we have the following proposition .
Note that we assume that the associated one-sided Riemann problem is uniquely solvable. Since the current paper is mainly concerned with a numerical algorithm for high order numerical boundary conditions, we leave aside for the moment the investigation of the rigorous mathematical theory.
2.4. One-sided Riemann solver and one-sided GRP solver
So-called solvers refer to the processes numerically solving the corresponding problems. Standard numerical Riemann solvers can be found in  and the generalized Riemann problem (GRP) solver in [3, 4]. The one-sided solvers proposed here are associated with the Riemann solver  and the GRP solver [3, 4]. These solvers aims (i) to provide all physical variables on the boundary ; (ii) to apply the inverse GRP to inspect the interaction of boundary and initial data.
Note that the boundary value that we obtain is not necessary to be continuous with the initial data at the origin
. If so, the solution is discontinuous. Such observation is heuristic when dealing with the interaction between shocks and solid boundaries. Besides, such a process provides several indications:
The compatibility of the resulting boundary data and the initial data determines the regularity of flows (solutions) around the origin locally. The one-sided Riemann solution is a key clue to the well-posedness. The one-sided Riemann solver aims to find numerically.
The one-sided GRP solution depends on the associated Riemann solution, and the corresponding GRP solver aims to find the value and helps to build high order numerical schemes.
2.5. One-sided GRP solver in two dimensions
We extend the one-sided GRP solver to two dimensions in this part. Suppose we have a boundary which is independent of time. For the boundary conditions that depend on the time such as a piston problem, we refer to  for the associated GRP solver. Our strategy includes the following steps: we first solve a normal one-sided Riemann problem at any fixed point on the boundary along the normal direction, namely,
where the boundary value is prescribed to be . Denote by the unit normal vector of . For the presentation simplicity, the boundary is set along the -axis, thanks to the Galilean invariance for fluid dynamical systems. Then (2.21) can be transformed to solving the following normal generalized Riemann problem along the -axis,
where , and are the Galilean transform of , and , respectively. After resolving , we transform back to obtain . The same as the 1-D case, we solve the normal conservation law at
to obtain the normal Riemann solution . Then we solve the following IBVP,
at , where the term is a fixed value with the instantaneous value obtained from (2.23) and interpolated from , reflecting the tangential effect along the boundary . Then the 2-D one-sided GRP solver follows exactly the same as the 2-D GRP solver, one can find more details in .
2.6. High order numerical boundary conditions
Once the one-sided GRP solver is available, the boundary data can be approximated with second order accuracy and the boundary volume can be dealt with as the ordinary control volume. That is, if at moment , and are known, then the boundary flux is approximated in a common way,
Furthermore the integral of source term can be evaluated using the interface method too,
Analogously, we deal with the 2-D case.
3. Application to Gas Dynamical Systems
In this section we discuss the one-sided Riemann problem (RP) and the one-sided generalized Riemann problem (GRP) for gas dynamical systems. We first discuss one-dimensional case in the form (2.1) with
and is a problem-dependent source term. In particular, for nozzle flows takes the form
where are the density, velocity and pressure of the fluids, respectively. is the cross-section area of the duct. is the total energy, the internal energy is given by the equation of state (EOS) . Note that (3.2) also includes the case of radially symmetric flows . The reactive Euler flows  can be treated similarly.
System (3.1) has three eigenvalues
The one-sided Riemann solver for Euler equations. As we pointed out in the last section, the one-sided Riemann problem plays a role in the justification of local well-posedness besides its numerical value. This problem is formulated as
The method solving this problem (3.4) follows the one for the classical Riemann problem. We fix the wave curve associated with from the state in the phase space, –space, and then investigate the solvability for the prescribed data . It is easily checked that the solvability of such a problem is up to the following two conditions:
There is an intersection point of and ;
The dimension .
The following are two typical examples, see Fig. 3.1.
Prescribed velocity . Then we find a point on so that is fixed. Certainly as , and so that one additional condition is needed. For instance, we can supplement the gas density as the gas property on the boundary . However, as , we need to check the Mach number . If , the boundary condition is not suitably prescribed.
Given upstream Mach number . The given value actually implies
We look for its intersection point with to find and then using the equation of state (EOS).
In summary, we can investigate the one-sided Riemann problem to identify that whether the upstream flow is supersonic or not as well as make clear the correct prescription of boundary conditions.
One-sided GRP solver. The one-sided GRP solver serves to solve (2.19) numerically. Assume that is (or approximated by) a smooth function with regular limiting values
Based on the corresponding one-sided Riemann problem, we can obtain the limiting value on the boundary . Essentially there are two versions in analogy with the standard GRP solver: An acoustic version and a nonlinear version.
Specified to the Euler equations, the one-sided GRP solver is implemented as follows.
Judge from the associated one-sided Riemann solution whether there emit strong waves in order to determine to use the acoustic or nonlinear GRP solver.
The acoustic GRP solver is the same as the linear case above.
The nonlinear GRP solver consists of two cases: a supersonic upstream flow and a subsonic upstream flow.
A supersonic upstream flow. All conditions are given at boundary .
A subsonic upstream flow. We apply the same procedure of the standard GRP solver  and naturally derive the one-sided relation
where the coefficients and are fully determined by the values and the slope value , the detailed expressions can be found in .
We are in position to compute the partial derivative values and from (3.7). If the boundary condition is given as and subsequently , then follows by the linear relation (3.7). As for the density derivative , we have
on the boundary from the EOS. Here is the local sound speed.
2-D one-sided GRP solver for Euler. The 2-D compressible Euler equations can be written as
where and represent the density, velocity, velocity, pressure and total energy, respectively. The 2-D one-sided GRP solver is the practical combination of the above 1-D one-sided GRP solver and the 2-D GRP solver [4, 12]. A key point is that the transversal effect is included in the solver development .
4. Implementation of the one-sided GRP scheme
4.1. Brief summary of the one-sided GRP scheme
So far, the one-sided GRP solver is developed to suit the GRP scheme near the boundary. The boundary control volume is then treated the same as the interior control volumes in the finite volume framework. We discretize the domain by equally computation mesh size and set . The cell represents the left boundary cell centered at and the cell represents the right boundary cell centered at , as shown in Fig. 4.1.
Since a standard finite volume method, such as the GRP method in , can be applied over computational cells in the interior domain, we only focus on the boundary cell . The one-sided GRP scheme at the boundary cell assumes the piecewise linear data
The vector is the constant slope of over cell at time with the time step size. To obtain the second order accuracy, the mid-point value is used
in the resolution of numerical flux and the source term discretization. We apply the 1-D one-sided GRP solver in the following steps.
Step 1. Given the piecewise linear initial data (4.1), approximate the mid-point value as follows,
The computation of is the main ingredient of the one-sided GRP scheme. The value is the local solution at to the following one-sided Riemann problem:
which can be solved by an exact or approximate one-sided Riemann solver . Here we apply the one-sided GRP procedure (3.7)-(3.9) to obtain the instantaneous value on the boundary , and then approximate using (4.3).
Step 2. Evaluate the next time values by using the following formula
where the source term is discretized with the mid-point rule in time and the trapezoidal rule in space.
Step 3. In order to suppress local oscillations as discontinuities are present near the boundary, we update the slope by using the following monotonicity algorithm limiter
In two-dimensional computations, we take rectangular meshes , , , as an example for simplicity, here centered at the grid point . The finite volume formula is applied over all cells ,
The initial data at time is expressed as bilinear functions
The values and can be analytically derived by the resolution of a local quasi 1-D GRP solver at each interface. The one-sided GRP solver is applied on the boundary. Then we can take the same procedure as that for 1-D case to implement the finite volume scheme.
4.2. Numerical Examples
We will present several numerical examples to validate the performance as the one-sided GRP solver is used. The examples include the interaction of shocks with solid boundaries, the radially symmetric flows, the nozzle flows, the Mach reflection of shock and the forward facing step problem.
Example 1. The scalar equation. We first use the Burgers equation to test the performance of the one-sided GRP solver. Consider the following initial-boundary value problem for the Burgers equation
The solution has an explicit formula
A compressible wave propagates to the left and forms a shock at on the boundary. As , a shock from propagates to the right. We compute the solution using the GRP scheme with the reflective boundary condition and the one-sided GRP solver, respectively. The solution is plotted from time to time in Fig. 4.3, from which it is observed that the one-sided GRP solver gives very sharp resolution of the singularity point , compared with the reflective boundary condition treatment.
Example 2. A single shock interaction with a solid boundary We test the example that a single shock wave interacts with a solid wall to verify the numerical performance of the one-sided GRP solver. The computational domain is [0,10] where the boundary is at . A left-propagating shock wave is initially positioned at . We take and the initial data is set to be
A reflected shock wave is observed when the output time is set to . We compare the results by two different boundary condition treatments: the traditional reflective boundary condition and the one-sided GRP solver. From Fig. 4.4, one can observe that the result obtained by the one-sided GRP solver is more stable and has less oscillations near the boundary. We further plot equally-distributed density contours from time to time at every time interval in Fig. 4.5, one can see again that the one-sided GRP solver gives very sharp resolution at the interaction point of the shock with the boundary.
Example 3. The Woodward-Colella problem. This is a classical interacting blast wave problem with the gas initially at rest and . The density is everywhere unit, the pressure is for and for , while it is only for . The solid-wall boundary conditions are prescribed at both ends. We compare the results of the reflective boundary condition treatment with that of the one-sided GRP solver. The CFL number is 0.6. The output time is set to . The numerical results for both boundary condition treatments are shown in Fig. 4.6 with 400 cells and 800 cells, respectively. It can be seen that the one-sided GRP solver is effective and robust for the blast wave problem.
Example 4. The nozzle flow. The nozzle flow problem is a classical quasi one-dimensional problem. Consider a flow in a converging-diverging nozzle occupying the domain . The cross-sectional area function of the duct is given by
with and . The governing equations are the Euler equations with geometric source term (3.1),(3.2). Set as the entrance of the duct and as the exit. We are concerned with the present boundary treatment to attain the steady state solution.
Two types of steady states are discussed: A continuous steady state and a discontinuous steady state containing a standing shock wave. The initial data for both cases can take as
where and are parameters to be determined, is a constant value determined by the steady solution at . In the previous study [3, 4], the inflow density, velocity and pressure are assigned to the inflow boundary condition, the outflow pressure is assigned as the outflow boundary condition. Here we apply the one-sided GRP solver to test its ability of attaining steady solutions.
For the first case, we set and in (4.13). This produces an isentropic continuous steady solutions which is defined by
in which the Mach number is determined by through the algebraic relation
In this case, the flow is transonic across the throat at the position . Thus the inflow boundary condition at the entrance should be prescribed by
While at the exit , the flow is supersonic and no boundary condition is needed. The computational result is given in Fig. 4.7 where 22 cells are used. The CFL number is 0.6 and the output time is . The solution obtained by implementing the one-sided GRP solver converges to the exact steady one and is comparable with the result obtained in .
For the other case, where the steady solution contains a standing shock wave, we set and in (4.13) to get the initial data. In this case, the flow jumps from supersonic to subsonic after passing the standing shock wave. As the outflow is subsonic in this case, both inflow boundary condition and outflow boundary condition should be imposed. The inflow boundary condition is at the entrance and the outflow boundary condition is at the exit . The computational result with 22 cells is given in Fig. 4.8. The CFL number is 0.6 and the output time is . The solution obtained by taking the one-sided GRP solver matches well with the exact solution.
Example 5. The spherical symmetric shock interaction problem. We test the one-sided GRP solver for the simulation of the spherical symmetric flows where a spherical shock wave interacts with the symmetric center. The initial data is taken to be