Stable, entropy-pressure compatible subsonic Riemann boundary condition for embedded DG compressible flow simulations

One approach to reduce the cost to simulate transitional compressible boundary layer flow is to adopt a near body reduced domain with boundary conditions enforced to be compatible with a computationally cheaper 3D RANS simulation. In such an approach it is desirable to enforce a consistent pressure distribution which is not typically the case when using the standard Riemann inflow boundary condition. We revisit the Riemann problem adopted in many DG based high fidelity formulations. Through analysis of the 1D linearised Euler equations it is demonstrated that maintaining entropy compatibility with the RANS simulation is important for a stable solution. The maintenance of Riemann invariant at outflow leaves one condition that can be imposed at the inflow. Therefore the entropy-pressure enforcement is the only stable boundary condition to enforce a known pressure distribution. We further demonstrate that all the entropy compatible inflow Riemann boundary conditions are stable providing the invariant compatible Riemann outflow boundary condition is also respected. Although the entropy-pressure compatible Riemann inflow boundary condition is stable from the 1D analysis, 2D tests highlight divergence in the inviscid problem and neutrally stable wiggles in the velocity fields in viscous simulations around the stagnation point. A 2D analysis about a non-uniform baseflow assumption provides insight into this stability issue (ill-posedness) and motivate the use of a mix of inflow boundary conditions in this region of the flow. As a validation we apply the proposed boundary conditions to a reduced domain of a wing section normal to the leading-edge of the CRM-NLF model taken out of a full 3D RANS simulation at Mach 0.86 and a Reynolds number of 8.5 million. The results show that the entropy-pressure compatible Riemann inflow boundary condition leads to a good agreement in pressure distribution.

enforce a consistent pressure distribution which is not typically the case when using the standard Riemann inflow boundary conditions. We therefore revisit the Riemann problem adopted in many DG based high fidelity formulations.
Through analysis of the one-dimensional linearised Euler equations it is demonstrated that maintaining entropy compatibility with the RANS simulation is important for a stable solution. It is also necessary to maintain the invariant for the Riemann outflow boundary condition in a subsonic flow leaving one condition that can be imposed at the inflow boundary. Therefore the entropy-pressure enforcement is the only stable boundary condition to enforce a known pressure distribution. We further demonstrate that all the entropy compatible inflow Riemann boundary conditions are stable providing the invariant compatible Riemann outflow boundary condition is also respected.
Although the entropy-pressure compatible Riemann inflow boundary condition is stable from the one-dimensional analysis, two-dimensional tests highlight [9] and convergence [10,11]. Compared with the strong enforcement, the flow quantity distributions are usually not directly set in the weak enforcement.
In addition, in the compressible flow regime, subsonic flows and supersonic  waves are stabilized by a favorable pressure gradient and destabilized by an adverse pressure gradient while both pressure gradients destabilize the crossflow waves [13,14]. Further the pressure load is usually well captured by the Euler or RANS simulation (at least for lift prediction) and the pressure distribution does not typically vary much over the boundary layer. The pressure distribution is 70 therefore considered a reliable quantity of interest from the lower fidelity models and is an important property to be maintained in the reduced domain.
To achieve the pressure compatibility in the reduced domain, the corresponding subsonic inflow boundary condition needs to be constructed. However, although this construction can follow several approaches, not all of them 75 are appropriate. The ideal boundary conditions enable the numerical simulation to converge in a stable way, inferring that the partial differential equation problem is well-posed and the numerical discretisation is stable with the adopted boundary conditions.
For the stability analysis, Giles [15] and Darmofal et al. [16] have developed 80 one-dimensional (1D) linear analysis (or eigenmode analysis) to evaluate decay rates of the disturbances. Since the decay rates depend on the adopted boundary conditions, the performance of boundary conditions are therefore obtained. This method can also be used to check the well-posedness of the boundary condition enforcement. Nordström et al. [12,17] used the energy method to analyze the 85 linearized system in a control volume. This method can be used to understand the number of conditions to make the problem well-posed as well as examine suitability of certain combination of boundary conditions. These works provide insightful analysis and set up a useful framework, which are wildly used to establish boundary conditions for compressible flow simulations. 90 As for the boundary condition stability analysis of DG approximations, it is less studied compared with those based on the conventional finite difference discretisation. Hindenlang et al. [18] examined the linear stability of the wall boundary condition for the inviscid flows using the sufficient condition developed by Gassner et al. [19]. This study focused on the wall boundary condition while 95 reasonably assumed the conditions on other boundaries are stable. However, since it is the combination of boundary conditions that determines the evolution of the system, the satisfaction of this condition on individual boundaries does not necessarily guarantee a non-singular system. To obtain a specific analysis within the conventional method, where the disturbance is introduced internally, 100 a uniform baseflow is sometimes needed to compute the integration over the boundaries of the control volume. This can limit its application to problems with large gradient, and therefore an extension of this methodology is required.
In this paper we therefore revisit the Riemann problem to see how best to enforce pressure compatibility at the subsonic inflow boundary. In Section 105 2, one-dimensional (1D) Euler equations are linearized in a DG element, and the systems for the inflow boundary conditions to be stable are derived from two perspectives. An entropy-pressure compatible inflow boundary condition is selected from the full range of possibilities in Section 3 . Section 5 extends the current analysis to multi-dimensions without assuming an uniform baseflow 110 assumption. This method is then used to understand the instability issue in adopting a stale 1D entropy-pressure inflow in the region of a stagnation point flow. Finally, a transonic flow in a reduced domain taken out of a full threedimensional (3D) RANS simulation is computed for validation.

115
To analyze the applicability of weak boundary conditions, we now consider the stability of the linearized Euler system subject to a discontinuous Galerkin (DG) discretisation. A common practice for DG approximations of the compressible flow equations is to impose boundary conditions using one-dimensional characteristic line theory. We therefore follow this practice and derive the linearized DG approximation for 1D Euler equations, whose conservative form where Q is the state vector of conservative variables, F is vector for inviscid flux, that is In the above ρ is the density, u is the velocity in the x-direction, p is the pressure and E is the total energy, which can be expressed under the assumption of a perfect gas law as where γ is the specific heat ratio.
In the DG approach, the approximation in each element are independent and boundary conditions are enforced weakly in each element in the form of a numerical flux. Therefore, to derive the smallest linear system for stability analysis which preserves all the necessary boundary conditions, we consider the domain as consisting of a single element Ω, as shown in Fig. 1. Eq. (1) is first multiplied by a test function φ and then integrated in the domain by part, which gives For simplicity we next consider a piecewise constant approximation in the Imposing this piecewise constant approximation to Eq. (3) allows us to analytically evaluate the first term and directly sets the third term to zero. Finally replacing the boundary flux F by a numerical flux denoted byf , Eq. (3) becomes where ∆x = x e − x w is the length of the 1D domain and the superscript w and e denote the west and east boundaries of the element. The numerical flux is typically computed through a Riemann solver, and takes the form where at the interface between two adjacent elementsQ int andQ ext represent the internal state and external states. For the numerical flux at the solution domain boundary, where there is no adjacent element, a ghost stateQ gh is introduced instead of the external state so that we havẽ In the standard use of a Riemann boundary condition the external or ghost state is independent of the interior state and the solution of the Riemann solver will provide a known boundary stateQ b , i.e.
In this circumstancesQ b is not known explicitly but rather inferred througĥ Since the numerical flux is constructed with the reference state and the internal state, the linearized system for stability analysis can be generated by introducing disturbances on both states, as is shown in Fig. 2. The linearized form of Eq. (5) (see Appendix A for full details), is given by where δQ int is the time-dependent disturbance on the internal state and δQ ref is the constant disturbance on the reference state. (The disturbance on the west boundary is assumed equals to that on the east boundary.) Although the independent variables in Eq. (9) can be in any form for the 1D analysis it is convenient to use the characteristic variables which are where c = γp/ρ is the speed of sound, S, which is the measure of entropy and takes the form By assuming the baseflow is steady and therefore uniform because of onedimensionality, the 1D linearized system in characteristic variables takes the where the coefficient matrices are defined as In the above, C 1 relates to the uniform baseflow quantities and can be further evaluated as: which is similar to the Jacobian matrix ∂f /∂Q, and therefore the eigenvalues are the well-known u + c, u, and u − c. The remaining coefficient matrices are which contain the normalized coefficients about how the boundary condition will 130 respond to the disturbance of the internal and reference states, respectively.
For the coefficient matrix C int , we assume it can be diagonalised so that The second condition is equivalent to a opposite relation for the second part of the coefficient matrices as If this condition is not satisfied, the solution will converge to This solution fails to track the reference value, and therefore is non well-posed.
Another possible scenario is when C int is singular so that some eigenvalue is zero, e.g. λ i = 0, the i-th component solution becomes  According to the above analysis, the two requirements to form a stable system are provided by Eq. (18), based on the eigenvalue analysis of the coefficient matrices. We therefore next consider how the coefficient matrices are constructed for specific Riemann boundary conditions.

Characteristic treatment for Riemann boundary conditions
where subscript l and r denote the left (ghost) state and right (internal) state, 155 respectively. The simplest boundary condition one can enforce is to specify the reference state to be a given ghost state independent of the internal values (as is in Fig.   2 which leads to the so called Riemann inflow [16]. We will also refer to this boundary condition as an entropy-invariant compatible inflow in the current work since the entropy (as well as the other incoming Riemann invariant) are directly enforced on the boundary state by an exact Riemann solver.

160
Since this inflow boundary condition only enforces entropy and the incoming Riemann invariant, there is no guarantee for pressure compatibility with the outer solution. To complete the stability analysis, the boundary condition at the outflow also needs to be provided. In what follows we will enforce the invariant compatible condition at the outflow in all of the following cases. This inflow-outflow combination, denoted by superscript SI, leads to the following coefficient matrices By substituting Eq. (24) and (25) into Eqs. (12) and (13), the coefficient matrices becomes where the second condition in Eq. (18) is directly satisfied. The eigenvalues of C SI int are −(c + u), −u, and −(c − u) while eigenvalues for C SI ref are the same but negative values. The physical interpretation of these eigenvalues is that when the internal state is disturbed, the disturbance can be decomposed onto the three waves, which leave the domain at the speeds c + u,u, and c − u. The minus signs represent the decaying feature of the disturbances. On the other hand, if the reference state is disturbed, the disturbance can also be projected onto the waves that come into the domain at the positive speed c + u, u, and c − u. Since the eigenvalues for C SI int are all negative, the first condition in Eq. (18) is also satisfied. Therefore, the solution to the linearized system in Eq. (11) is stable 170 at the steady state using these boundary conditions.
We note from inspecting Eqs. (24) and (25) that the following relations hold and The stability of this system can be analyzed directly through the eigenvalues of 180 C int . Alternatively, we can also examine the stability by constructing C ref and using the opposite relation.
A special case that is worth further discussion is when one or more eigenval- Similar to the analysis for strong boundary conditions by Darmofal et al.

200
[16], the coefficient matrix C 2,int contains the reflection properties of the adopted boundary conditions. However, here, it is the component matrices ∂Û w b /∂Û int and ∂Û e b /∂Û int rather than C 2,int should be checked to avoid entry cancellations. The zeros of ∂Û w b /∂Û int Therefore, the combination of entropy-invariant compatible inflow and invariant compatible outflow is fully non-reflective in the normal direction.

Stability analysis for pressure compatible inflows
We now analyzes the eigenvalues for C int to determine a stable pressure 210 compatible inflow for the DG solver. Three possible candidates are constructed and their stability performances analyzed. We only examine the coefficient matrix C int since the opposite relation in Eq. (19) always hold.

Entropy-pressure compatible inflow
Neglecting the superscript w, the entropy-pressure compatible inflow at the west boundary in Fig. 1 is given by which directly implies The third condition which needs to be specified according to the characteristic relations in Eq. (22) Then the boundary state velocity and right-propagating characteristic are computed as Therefore, the boundary state Riemann variables take the form As mentioned earlier, using the invariant compatible Riemann outflow, the with which the eigenvalues for C SP int are Since all of the eigenvalues are negative for subsonic flows, the combination of 215 the entropy-pressure inflow and the invariable compatible outflow leads to a stable system.
It is worth noting that the entropy-pressure compatible inflow is able to enforce more conditions in multiple-dimensional cases, where the tangential velocity components are enforced by the reference state. The density compatibility 220 in Eq. (33) makes the boundary state tangential momentum components also match the reference state.

Velocity-pressure compatible inflow
Following the same notation, the velocity-pressure compatible inflow satisfies Then the boundary state density is computed using the left-pointing characteristic as Therefore, the boundary state characteristic variables are given by which leads to and the eigenvalues of C U P int are where λ U P 2,3 are purely imaginary for subsonic flows. The first order linear analysis, indicates that part of the internal disturbance will keep oscillating and 225 therefore the velocity-pressure compatible inflow is not desirable for a steady state solution.

Momentum-pressure compatible inflow
Similarly, a momentum and velocity boundary state can be enforced by the conditions By introducing a auxiliary variable, η, the boundary state density, velocity, and speed of sound take the form Therefore, the boundary state characteristic variables are The second part of coefficient matrix is therefore given by Subsequently the eigenvalues of C M P int are where Since the second and third eigenvalues have complicated expressions, their values are examined numerically at different subsonic Mach numbers M a, as 230 shown in Fig. 4. We observe that these two eigenvalues are always real and positive, leading to a unstable solution and so the adopted boundary conditions are unstable.

Stability analysis of inflow entropy modification
In Section 3 we demonstrated that the only stable inflow boundary condi- imply the entropy is a essential quantity, and violating entropy compatibility is likely to cause problems. We therefore provide more detailed analysis for 245 entropy modification at inflow.
We continue to follow the setting in Fig. 1 where the west boundary is the inflow, and again the invariant compatible outflow at the east boundary is adopted. Since only the outward-propagating Riemann invariant may be used to construct the boundary state at the inflow, the second part coefficient matrix takes the form where a and b are entries representing the contribution of the internal state.
Then the characteristic polynomial for the coefficient matrix is where λ 1 = −(c + u) is always negative and λ 2,3 normalized by speed of sound       on the right wall. We next construct a 2D linearized system using the DG approximation of the 2D Euler equations, which take the conservation form where Q is the state vector of 2D conservative variables, F 1 and F 2 are the vectors of the inviscid flux In the above, u 1 and u 2 are the velocity components in the x 1 and x 2 -directions, respectively, and E = 1 γ−1 p + 1 2 ρ(u 1 2 + u 2 2 ). Others symbols are the same as in Eq. (1).
Following a similar analysis as performed in Section 2 and Appendix A, the domain Ω is approximated by a single quadrilateral DG element as The expanded form of Eq. (59) on the squared domain in Fig. 6(a) is where the subscript n and s denote north and south, respectively. Due to the presence of the wall the baseflow is no longer uniform so that a two-dimensional approximation of the following form is considered where P 1 and P 2 are expansion order in the orthogonal directions, and typically are set to be equal to each other P 1 = P 2 = P . In the following analysis we choose to adopt a Lagrange nodal basis function for both φ p and φ q of the form , otherwise, where x p represents the Gauss-Lobatto-Legendre points including both endpoints of the interval. For clarity, we use the semi-discrete form by only numer-ically integrating the boundary terms where the numerical flux going through the wall,f e 1,i , is constructed using the internal state and no-penetration condition provided through the reference statê where the contributions by internal state disturbances and reference state disturbances are respectively included in the right-hand side terms RHS int and RHS ref , are given by and In the above, due to the presence of the wall, neither the baseflow nor the disturbances have a uniform distribution. Therefore the right-hand side terms cannot be further simplified to form two coefficient matrices as we have derived in the 1D analysis. So from the conventional point of view, the eigenvalue analysis for the coefficient matrix only with respect to the internal disturbance cannot be applied to this system. However, the instability issue can be shown by checking the responses to the reference state disturbances through the derivative matrix, i.e. ∂V l b,i /∂V l ref,i (l = w, s, n) which arises in each term of RHS ref . Considering the boundary conditions adopted, the boundary states for the east, north, and south boundaries are given bŷ where the subscript "n/s" denotes either north or south boundary, and the velocity tangential to the boundary is enforced to be the upstream value according to the feature of split Riemann problem [21]. The respective derivative matrices of boundary state with respect to the internal and reference states are In this case the non-zero entry in the second column enables the disturbance on the reference state u 1 -velocity component go into the domain to determine the 310 field and therefore the problem becomes well-posed.

Results
Having identified out the stable entropy-pressure compatible inflow and understood its potential singularity issue around a stagnation point, we next provide an example for pressure compatible simulation in the reduced domain on 315 the CRM-NLF model [22]. The Navier-Stokes equations are solved using the open-source spectral/hp framework Nektar++ [23,24].
To generate the background fields, a transonic RANS simulation over the full wing-fuselage geometry is first carried out under the conditions provided in Table. 2, where α is the angle of attack, L is the reference length. Fig. 7 shows the pressure distribution on the surface where two shocks on the wing can be observed. We therefore set the outflow boundary of the reduced domain at the upstream of the shock. The reduced domain is obtained on the slice in the wall-normal direction passing the leading-edge of section D of the CRM-NLF model [22].  In Fig. 9 we compare the pressure distributions when using different inflow boundary conditions. The result for entropy-pressure compatible inflow is given in Fig. 9(b) where the contour lines perfectly match the pressure contours of 340 the RANS field. Fig. 9(c) and Fig. 9(d) provide the pressure distributions without pressure compatibility enforced. Significant mismatching in contour lines is observed in the result for the entropy-invariant compatible inflow in Fig. 9(c), where the invariant compatible outflow rather than pressure outflow needs to be adopted for a stable simulation because of the large pressure incompatibility.  periodic boundary condition is adopted in the x 3 -direction for the 3D simulation. Fig. 9(d) shows the pressure distribution for using the entropy-enthalpy compatible inflow. This result is better than that for entropy-invariant compatible inflow but some contour lines still unable to match the RANS data.

Conclusion
To carry out high fidelity simulation in a reduced domain which is embedded in a outer domain where a lower fidelity model such as RANS is applied, appropriate boundary condition enforcement is required. This is particularly the case when the high fidelity simulation is compressible but subsonic, since not all the data interpolated from the outer simulation can be enforced at the inflow bound-365 ary of the reduced domain. Since the pressure load is typically well predicted by the low fidelity model and does not vary much over the boundary layer, the pressure is considered a reliable quantity to maintain in the reduced domain.
We have therefore determined the stable enforcement to achieve pressure compatibility at subsonic inflow boundary for DG compressible flow simulations.  Since theQ e ref,i will not be influenced by the outer simulation, we have δQ e ref,i = 0. The consequent linearized system for Fig. 6(a) is derived as In the above RHS ref has two less terms than RHS int .