1. Introduction
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[27], mathematical justifications (wellposedness arguments) [16] and numerical treatment. Physical considerations prescribe boundary data for a part of physical variables based on specific problems such as the solidwall boundary condition; mathematical wellposedness 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 [19]. In [27], characteristic method is used but restricted to first order accuracy for smooth flows. Although the resulting schemes may be wellimplementable, 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 [27]. 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 wellposed. 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].We associate this issue with the socalled onesided generalized Riemann problem (GRP). As motivation, we take a look at the initialboundary value problem for the Burgers equation [2, 1],
(1.1) 
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
(1.2) 
, . For this case, we have a shock solution,
(1.3) 
, . 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 [2], 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 [6] and the solid body floating in the air [13]. 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 onesided GRP. It is an initialboundary value problem rather than a purely initial value problem. Correspondingly, a numerical method to solve this problem is called a onesided GRP solver. Such a study has twofold 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 onesided GRP formulated here and they have been refined to put into the simulation of fluid flows with moving boundaries [14, 13]. In [12], the similar idea was employed if the flow is smooth and consistent with the inverse LaxWendroff method [29]. Note that since the twostage 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 onesided GRP problem. In Section 3 we discuss numerical boundary conditions for compressible Euler equations via the passage of onesided 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 onesided GRP solver for hyperbolic balance laws
In this section we formulate onesided 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 initialboundary value problem for compressible fluid flows in [25] 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 initialboundary value problem is formulated below as the onesided generalized Riemann problem (OSGRP).
Consider hyperbolic balance laws
(2.1) 
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 [10]. It is assumed to be hyperbolic in the sense that the Jacobian of has
real eigenvalues
with a complete set of associated eigenvectors
,(2.2) 
Each is genuinely nonlinear or linearly degenerate in the sense of Lax [20].
We focus on the left boundary . The right boundary is treated similarly. We emphasize that the free boundary problem can be studied too [13]. On the boundary , the data is imposed as
(2.3) 
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
(2.4) 
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 ,
(2.5) 
Then we propose the following assumption.
Assumption. The problem (2.1) –(2.3) is wellposed at least locally so that
(2.6) 
in some “appropriate” sense.
This assumption is very “rough” and understood in certain intuitive way. A primary judgement of the wellposedness boils down to the solvability of the following onesided Riemann problem, while the dynamics is dependent on the onesided 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 onesided 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
(2.7) 
The characteristic decomposition tells that
(2.8) 
where , is a lefteigenvector associated with the eigenvalue . The solution formula is
(2.9) 
To obtain the solution values on the boundary , we have from (2.8)
(2.10) 
Hence the boundary value of can be obtained by solving the following system
(2.11) 
Indeed, the wellposedness of (2.7) depends on the solvability of (2.11). That is,
(2.12) 
where ,
, 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 onesided 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
(2.13) 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
(2.14) 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 onesided Riemann problem and onesided generalized Riemann problem, respectively.
2.2. Onesided Riemann problem
The onesided Riemann problem is motivated from the piston problem [8] and formulated in [13]. Here we formulate this problem for hyperbolic conservation laws
(2.15) 
The boundary data is prescribed as
(2.16) 
for some , where the operator prescribes certain physically meaningful values to partial state components. Corresponding to (2.1)(2.3), and .
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 onesided 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 wellposed. Hence this onesided Riemann problem plays a role in checking whether the boundary conditions are correctly prescribed.
Another role of the onesided 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,
(2.17) 
Then it is necessary to check whether there are exactly characteristics leaving from the boundary , similar to the linear case.
(2.18) 
That is, the dimension of manifold Just like the case for the Burgers equation, this is not necessary true. Hence the solvability of onesided Riemann problem is a necessary to judge the wellposedness of initialboundary value problem for (2.4).
2.3. Onesided 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 onesided 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 onesided generalized Riemann problem (GRP) is reformulated here as the initial and boundary value problem,
(2.19) 
where is smooth, and is measurable. This is associated with the onesided Riemann problem above. Similar to the interrelation between the standard generalized Riemann problem and the associated Riemann problem, we have the following proposition [4].
Proposition 2.1.
Note that we assume that the associated onesided 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. Onesided Riemann solver and onesided GRP solver
Socalled solvers refer to the processes numerically solving the corresponding problems. Standard numerical Riemann solvers can be found in [31] and the generalized Riemann problem (GRP) solver in [3, 4]. The onesided solvers proposed here are associated with the Riemann solver [31] 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 onesided Riemann solution is a key clue to the wellposedness. The onesided Riemann solver aims to find numerically.

The onesided 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. Onesided GRP solver in two dimensions
We extend the onesided 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 [13] for the associated GRP solver. Our strategy includes the following steps: we first solve a normal onesided Riemann problem at any fixed point on the boundary along the normal direction, namely,
(2.21) 
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,
(2.22) 
where , and are the Galilean transform of , and , respectively. After resolving , we transform back to obtain . The same as the 1D case, we solve the normal conservation law at
(2.23) 
to obtain the normal Riemann solution . Then we solve the following IBVP,
(2.24) 
to obtain
(2.25) 
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 [24]. Then the 2D onesided GRP solver follows exactly the same as the 2D GRP solver, one can find more details in [24].
2.6. High order numerical boundary conditions
Once the onesided 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,
(2.26) 
Furthermore the integral of source term can be evaluated using the interface method too,
(2.27)  
Analogously, we deal with the 2D case.
3. Application to Gas Dynamical Systems
In this section we discuss the onesided Riemann problem (RP) and the onesided generalized Riemann problem (GRP) for gas dynamical systems. We first discuss onedimensional case in the form (2.1) with
(3.1) 
and is a problemdependent source term. In particular, for nozzle flows takes the form
(3.2) 
where are the density, velocity and pressure of the fluids, respectively. is the crosssection 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 [22]. The reactive Euler flows [30] can be treated similarly.
System (3.1) has three eigenvalues
(3.3) 
where is the local sound speed. All other properties of this system can be found in any textbook about gas dynamics, e.g. [8, 9, 3, 31].
The onesided Riemann solver for Euler equations. As we pointed out in the last section, the onesided Riemann problem plays a role in the justification of local wellposedness besides its numerical value. This problem is formulated as
(3.4) 
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
(3.5) We look for its intersection point with to find and then using the equation of state (EOS).
In summary, we can investigate the onesided Riemann problem to identify that whether the upstream flow is supersonic or not as well as make clear the correct prescription of boundary conditions.
Onesided GRP solver. The onesided GRP solver serves to solve (2.19) numerically. Assume that is (or approximated by) a smooth function with regular limiting values
(3.6) 
Based on the corresponding onesided 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 onesided GRP solver is implemented as follows.

Judge from the associated onesided 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 [5] and naturally derive the onesided relation
(3.7) where the coefficients and are fully determined by the values and the slope value , the detailed expressions can be found in [4].

2D onesided GRP solver for Euler. The 2D compressible Euler equations can be written as
(3.10) 
where and represent the density, velocity, velocity, pressure and total energy, respectively. The 2D onesided GRP solver is the practical combination of the above 1D onesided GRP solver and the 2D GRP solver [4, 12]. A key point is that the transversal effect is included in the solver development [21].
4. Implementation of the onesided GRP scheme
4.1. Brief summary of the onesided GRP scheme
So far, the onesided 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 [5], can be applied over computational cells in the interior domain, we only focus on the boundary cell . The onesided GRP scheme at the boundary cell assumes the piecewise linear data
(4.1) 
The vector is the constant slope of over cell at time with the time step size. To obtain the second order accuracy, the midpoint value is used
(4.2) 
in the resolution of numerical flux and the source term discretization. We apply the 1D onesided GRP solver in the following steps.
Step 1. Given the piecewise linear initial data (4.1), approximate the midpoint value as follows,
(4.3) 
The computation of is the main ingredient of the onesided GRP scheme. The value is the local solution at to the following onesided Riemann problem:
(4.4) 
which can be solved by an exact or approximate onesided Riemann solver [31]. Here we apply the onesided 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
(4.5) 
where the source term is discretized with the midpoint 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
(4.6) 
More details about the minmod function can be found in [5, 31].
In twodimensional 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 ,
(4.7) 
The initial data at time is expressed as bilinear functions
(4.8) 
The values and can be analytically derived by the resolution of a local quasi 1D GRP solver at each interface. The onesided GRP solver is applied on the boundary. Then we can take the same procedure as that for 1D case to implement the finite volume scheme.
4.2. Numerical Examples
We will present several numerical examples to validate the performance as the onesided 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 onesided GRP solver. Consider the following initialboundary value problem for the Burgers equation
(4.9) 
The solution has an explicit formula
(4.10) 
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 onesided GRP solver, respectively. The solution is plotted from time to time in Fig. 4.3, from which it is observed that the onesided 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 onesided GRP solver. The computational domain is [0,10] where the boundary is at . A leftpropagating shock wave is initially positioned at . We take and the initial data is set to be
(4.11) 
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 onesided GRP solver. From Fig. 4.4, one can observe that the result obtained by the onesided GRP solver is more stable and has less oscillations near the boundary. We further plot equallydistributed density contours from time to time at every time interval in Fig. 4.5, one can see again that the onesided GRP solver gives very sharp resolution at the interaction point of the shock with the boundary.
Example 3. The WoodwardColella 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 solidwall boundary conditions are prescribed at both ends. We compare the results of the reflective boundary condition treatment with that of the onesided 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 onesided GRP solver is effective and robust for the blast wave problem.
Example 4. The nozzle flow. The nozzle flow problem is a classical quasi onedimensional problem. Consider a flow in a convergingdiverging nozzle occupying the domain . The crosssectional area function of the duct is given by
(4.12) 
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
(4.13) 
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 onesided 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
(4.14) 
in which the Mach number is determined by through the algebraic relation
(4.15) 
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
(4.16) 
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 onesided GRP solver converges to the exact steady one and is comparable with the result obtained in [4].
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 onesided GRP solver matches well with the exact solution.
Example 5. The spherical symmetric shock interaction problem. We test the onesided 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
Comments
There are no comments yet.