Shock equations and jump conditions for the 2D Adjoint Euler equations

This paper considers the formulation of the adjoint problem in two dimensions when there are shocks in the flow solution. For typical cost functions, the adjoint variables are continuous at shocks, where they have to obey an internal boundary condition, but their derivatives may be discontinuous. The derivation of the adjoint shock equations is reviewed and detailed predictions for the behavior of the gradients of the adjoint variables at shocks are obtained as jump conditions for the normal adjoint gradients in terms of the tangent gradients. Several numerical computations on a very fine mesh are used to illustrate the behavior of numerical adjoint solutions at shocks.


Introduction
This work deals with some technical details of inviscid adjoint equations in the presence of shocks.Understanding how inviscid adjoint solutions behave at shocks is relevant for transonic and supersonic design and control applications [1], with potential applications ranging from drag or sonic boom reduction to fluidstructure and flutter control, but also for error analysis and correction [2].To understand the adjoint approach for shocked flows, three major points need to be addressed: 1.The correct formulation of the inviscid adjoint equations at shocks, which has to consider linearized perturbations to the shock location.This analysis was carried out in [3] [4] for the quasi-1D adjoint Euler equations and in [5] [6] for the 2D adjoint Euler equations.In both cases, it was shown that, for typical cost functions, the adjoint variables are continuous at shocks where they have to obey an adjoint boundary condition.The correct formulation and approximation of adjoint equations in flows with shocks has also been addressed in [7] [8] [9].2. The accuracy of discretized adjoint approximations with shocks.For numerical computation, the adjoint shock boundary condition is usually not explicitly enforced, but the discretized adjoint equations (either continuous or discrete) yield the correct solution as long as adequate levels of smoothing are applied, in such a way that the shock is spread over an increasing number of grid points as the mesh spacing is reduced [10] [11, 12]. 3. Assess the impact of the shock equations on practical applications such as aerodynamic design or error estimation, among others.In numerical implementation, especially in 2D, such shock conditions and increased smoothing are usually ignored under the assumption that they have little impact on both the adjoint solution and the sensitivities.However, results in [6] show that explicitly accounting for the adjoint shock condition can have a significant impact on both the sensitivities and the optimization procedure.For adjoint-based error correction with shocks, it was shown in [2] with 1D examples that the key to obtaining meaningful results is the use of discrete solutions with a well-resolved viscous shock.
The focus of the present work will be on point 1.While in quasi-1D the behavior of the adjoint variables and their derivatives can be completely established from the shock and adjoint equations [4] [13], a similar analysis for the 2D case is not yet fully available.It was argued in [5] that the adjoint variables, relative to an output function consisting of the integral of a function of the pressure along the airfoil surface, are continuous across shocks and that an adjoint boundary condition is required along the length of the shock.It was also claimed that gradients of adjoint variables are discontinuous, but no proof or conclusive numerical result was offered.A complete analysis of the 2D case involving was presented in [6], where it was shown that, for typical costs functions (aerodynamic lift and drag, for example), the adjoint variables are continuous across shocks and that the adjoint variables obey a differential equation along the shock.The former statement is cost function dependent.For example, entropy variables are adjoint states relative to a cost function measuring the net entropy flux across the domain boundaries (including the shock surface) [14] and are discontinuous across shocks.
Concerning the behavior of the gradients of adjoint variables at shocks, no further analysis was conducted in [6].A first step in that direction was taken in [13], which presented the jump conditions for the gradients of the adjoint variables at shocks.It was not possible to derive closed form results for all the derivatives, not even for normal shocks, but it was claimed, based on preliminary analysis and numerical computations, that the results were compatible with the derivatives of the adjoint variables along the shock normal direction being continuous and mostly vanishing, at least for normal shocks.However, normal shocks are difficult to realize in practice and none of the shocks presented in [13] were normal, so the conclusions offered were somewhat misleading.In [15], a series of detailed numerical experiments on very fine meshes confirmed that adjoint derivatives are indeed generally discontinuous across generic (oblique) shocks.
Here, we would like to complement the results of [13] [15] in two directions.First, after reviewing the derivation of the adjoint shock equations in the 2D case, we will expand the analysis of the jump conditions of the adjoint gradients to put it on firmer grounds.Second, we will examine several numerical solutions containing oblique and normal shocks to illustrate the differences between both cases in the light of the derived shock conditions.

Adjoint Equations for Shocked Flows
This section contains a review of some known technical details of 2D adjoint equations in the presence of shocks originally derived in [6] (see also [13]).We consider steady transonic inviscid flow past a body of surface S such that the flow contains a shock attached to S with profile  extending from b x (shock foot) to end x (shock tip).We define the following cost function where p is the pressure and S n is the unit normal vector to S (see Figure 1).When  () h p p , eq. ( 1) corresponds to the force exerted by the fluid on S measured along a direction d .The flow is governed by the nonlinear (steady) flow equations and the shock equations (the Rankine-Hugoniot conditions), which for a steady shock can be written as In eqs.( 2) and (3), are the velocity components in the local shock frame [16] (see Figure 1).For normal shocks,  0 t v on either side and thus the 3 rd equation in eq. ( 3) is trivially verified.
Let us now consider the problem of computing the sensitivity derivatives of J(S) given by eq. ( 1) with respect to, say, changes in the shape of S. A deformation of S along the local normal direction The adjoint approach for this problem uses two adjoint fields: a bulk adjoint       Linearizing eq. ( 6) using eqs.( 4) and ( 5) yields [6]  x is the jump across the shock at the shock foot.In eq. ( 7), the third term on the first line includes the effect of a linearized displacement  x n n S n t in the shock foot location (see Figure 3 and the appendix of [6] for a detailed derivation of this term).Integrating by parts the domain term in eq. ( 7) yields where the identity  ) , has been used to rewrite the wall integral, and  S denotes the far-field boundary.Inserting ( 8) into (7) and rearranging yields where the linearized wall boundary condition  has been used to rewrite the third term in (9).The final step is to integrate by parts (   is the local curvature of the shock profile) and Hence, we can write (10) as where we have used that on both sides of the shock.Gathering ( 9) and ( 11) and rearranging yields The whole point of the adjoint approach is to define the adjoint variables in such a way that  () JS can be computed independently of U and  .This can be achieved if the bulk adjoint state obeys the adjoint equation , with the following wall and far-field boundary conditions Along the shock, the adjoint variables obey the equations Hence,  is continuous across the shock and, in consequence, so is the tangent derivative   tg .The second equation in (15) then becomes an ordinary differential equation for the bulk adjoint  along the shock, where . Finally, the last equation in (15) gives an initial condition for eq.( 17) at the

Detached shocks
For sufficiently high incoming Mach number, transonic flow past blunt bodies contains detached (bow) shocks ahead of the body.The shock is normal directly in front of the body, and extends around it as a curved oblique shock.At a sufficient distance from the body the shock reduces to a Mach wave at points a x and b x .
In this case, the adjoint problem is essentially unchanged.Along the bow shock the adjoint state is continuous and obeys the differential equation ( 17), but the matching condition ( 18) is now missing, being replaced by a boundary condition that is trivially satisfied.

Extension to 3D
The extension to three dimensions poses no significant complications.In the shocked case, the shock is now a surface which, for the sake of the analysis, will be taken to be a single sheet attached to the wing surface along a curve  S and with an additional boundary curve  o along which the shock merges with the remaining smooth flow.
It can be shown that the adjoint state is still continuous across the shock, where it obeys the following shock equations on the shock surface, and along the shock foot  S .Here    S is the jump across the shock at the shock foot,   tg is the tangent gradient (the covariant derivative on the shock surface) and  ˆS n is the normal vector to the shock boundary curve  S but otherwise tangent to the shock surface.

Jump conditions for the adjoint gradient across a shock
The values of the adjoint variables and their derivatives at a shock are constrained by the shock equations ( 17) and ( 18), as well as the adjoint equations on either side of the shock, which we write in shock-oriented coordinates [16] as where F n and the superscripts ± are used to distinguish upstream and downstream quantities.Taking differences across the shock yields These equations, along with (17), were analyzed and tested on a numerical solution in [13], but the meshes employed were relatively coarse.A much thorough numerical assessment on very fine meshes was carried out later on in [15], where the following equations that follow from ( 22) by elementary manipulations, were shown to be obeyed to a high degree of accuracy.
Here, we would like to take a step forward and try to constrain as much as possible the values of the jumps of the adjoint derivatives across the shock.The resulting equations are exact, and numerical testing should be seen more as a test on the extent to which numerical solutions obey them than to a test of the equations themselves.It was argued above that the determinant   det( ) 0 n A at generic points of 2D shocks, so eq.( 21) can be used to solve for the normal adjoint derivatives in terms of the tangent adjoint derivatives Taking differences across the shock in (24) gives equations for the jumps of the normal derivatives that, using the RH conditions and the adjoint shock equation ( 17), can be written as The manipulations required to tackle (24)-( 26) are carried out with the aid of a symbolic manipulation software.From (25) it is possible to obtain two additional properties for the jumps Notice that, according to (25), the jumps in the normal adjoint derivatives are generally not zero, but they are zero if the adjoint solution is constant along the shock.This is in agreement with the fact that, for example, there is no discontinuity of the adjoint variables or their gradients across the fish-tail shocks in supersonic cases, since the adjoint variables are constant (actually zero) on either side of the shock (see section 4.1).
Finally, for costs functions that depend only on pressure (lift and drag, for example) we have   14 H [5] and, thus, using (17) we get (where it has been assumed here that the flow is homenthalpic) and as well as the detached shock in front of a sufficiently wide wedge.Normal shocks are also present in most supersonic inlets.On the other hand, the shock forming over an airfoil in transonic flow is generally not normal.
If flow separation after the shock is not considered, such shocks are then normal at their feet, but staying normal to some distance above the surface (a so-called normal shock with normal extension [17]) is only possible for one very specific upstream wall Mach number, which is equal to   1.662 For a normal shock, the shock equation ( 17) reduces to: according to our conventions in Figure 1) We can actually derive stronger conditions for the normal derivatives themselves by going back to the adjoint equations ( 21), setting  0 t v and solving for the normal adjoint derivatives in terms of the tangent adjoint derivatives using also eq. ( 31), which yields

Numerical Tests
We will now present several test cases in order to illustrate the behavior of typical numerical adjoint solutions at shocks.We do not attempt to investigate the conditions under which the numerical solutions converges to the analytic solution nor to compare different numerical schemes as regards to the behavior of the computed adjoint solutions, but rather to show how a typical, finite-volume adjoint solver behaves at shocks without enforcing any shock condition.
We will be using the unstructured, open source SU2 code [18] for numerical testing.The computations are carried out on a very fine unstructured mesh obtained from the basic triangular Euler mesh (Figure 4) available at the SU2 test case repository [19] by 5 rounds of uniform refinement.At each round, every edge is bisected and the resulting nodes are joined to form new triangles.In order to preserve the surface of the airfoil, a Bézierspline surface reconstruction on the basis of the previous mesh has been performed at each stage.The final mesh has 6400 nodes on the airfoil profile and 5.2×10 6 nodes and 10.4×10 6 triangular elements throughout the flowfield.The near wall distance is around 10 -5 chord lengths, which should be more than adequate to resolve the Euler flow and adjoint fields.Notice that this compares well to the current stat of the art on these type of mesh-converged Euler computations [15] [20].The flow and drag-based adjoint solutions are computed with the SU2 direct and continuous adjoint solvers using a central scheme with JST artificial dissipation.

Supersonic Flow
Our first test case is supersonic flow past a NACA0012 airfoil.Flow conditions are M = 1.5 and angle of attack α = 0°, which result in a detached bow shock ahead of the leading edge and two inclined fish tail shocks emanating from the sharp trailing edge (Figure 5).We now examine the behavior of the adjoint derivatives along a line crossing the bow shock as indicated in Figure 5 and Figure 7.The line is perpendicular to the shock, which is oblique at this location.The tip of the bow shock is a normal shock, which could be used to test the adjoint properties across such shocks, but the adjoint solution is contaminated by the presence of the adjoint singularity along the incoming streamline stagnation (the supersonic adjoint solution shows singularities along all three characteristics impinging the tip of the bow shock).
Figure 9, which plots the adjoint normal derivatives along the line (left panel).The largest jump is found for  3 , and the relative sizes of the jumps are in agreement with Eq. (29).
Figure 9 also checks on the right panel some of the properties of the jumps of the adjoint gradients presented in Eq. ( 26) and ( 27), as well as the adjoint shock equation (17).Bear in mind that what is being plotted is actually The plotted functions (36) are continuous across the shock, which indicates that the corresponding jumps are zero as required by Eq. ( 17), ( 26) and (27).

Transonic Flow
The second test case correspond to transonic flow past a NACA0012 airfoil with flow conditions M = 0.8, α = 1.25°.Under these conditions, the flow has a fairly strong shock on the upper side at x/c ≈ 0.64 and a weaker shock on the lower side (Figure 10).Contour plots of the adjoint variables are shown in Figure 11 and Figure 12.Notice that the adjoint solution is singular along the incoming stagnation streamline [5] [21], the wall [21] and also along a delta-like structure formed by the supersonic characteristic emanating from the shock foot and its reflection off the sonic line [22], but it is continuous across the shock and the sonic line.A zoom of the adjoint solution near the upper shock (Figure 12) shows that the adjoint derivatives are actually discontinuous (notice the abrupt change of direction of the adjoint contour lines of the y-momentum adjoint variable  3 ).We now examine the behavior of the adjoint derivatives along a line crossing the upper shock at right angles as indicated in Figure 10 and Figure 12.As explained above, this shock is not normal (the tangent velocity, though small, is definitely not zero, as can be seen in Figure 13).Figure 13 shows that the adjoint variables are clearly continuous across the shock, but their gradients are not (see the kink in the y-momentum adjoint variable  3 at the location of the shock).The discontinuity in the normal adjoint gradients is most clearly seen in Figure 14, which plots the adjoint normal derivatives along the line (left panel).The largest jump is found for  3 , while the jumps in the other normal derivatives cannot be appreciated at this level of resolution (and could very well be nearly zero; for  2 this is reasonable since, from Eq. ( 29),      2 [ ] 0 x n t , while for  1 and  4 it could be related to the fact that Hv and t v is relatively small).The right panel of Figure 14 checks the properties of the jumps of the adjoint gradients, Eq. ( 17), ( 26) and ( 27).The plotted functions are continuous across the shock, which indicates that the corresponding jumps are zero.

Normal shock with normal extension
For the angle of attack used in the previous case the upstream wall Mach number

Conclusions
In 2D, the adjoint variables of typical cost functions are continuous at shocks and obey a differential equation along the length of the shock.In this work, we have focused on determining the behavior of the gradients of the adjoint variables across shocks.It had been anticipated in [5], by studying the behavior of the Green's functions of the linearized Euler equations, that there is a discontinuity in the gradient of the adjoint variables at the location of the shock which was clearly visible in the numerical results of [15].
In order to analyze the behavior of the adjoint gradients at the shock, it proves convenient to use a frame of reference locally aligned with the shock curve, relative to which it is possible to decompose the adjoint gradients on either side of the shock into (locally) normal and tangent components.Since the adjoint variables are continuous at the shock, the tangent adjoint gradients are continuous as well, but the normal components need not.Combining the shock equation and the adjoint equations, which hold on either side of the shock, it is possible to obtain detailed predictions for the behavior of the gradients of the adjoint variables across the shock.These equations relate the discontinuity (jump) of the normal component of the adjoint gradients to their tangent component, but do not fix completely the value of the jumps.In fact, the equations allow continuous normal gradients provided that the tangent gradients vanish, as occurs for example across the fishtail shocks in supersonic flows.In 2D, the equations imply that, for normal shocks, 3 linear combinations of normal adjoint gradients are identically zero at the shock, while a fourth combination has a jump proportional to the discontinuity in the inverse normal component of the velocity Several numerical tests have been carried out to get a flavor of how numerical adjoint solutions behave at shocks.The computed solutions are continuous at the shock and present discontinuities in the gradient of the adjoint variables at the location of the shock.The solutions follow the shock relations to a reasonable degree, in agreement with the results presented in [15].
A separate question that we have not addressed here is that of the consistency of the adjoint solution for shocked flows in the limit of infinite grid resolution.In numerical computations, the internal boundary condition (the differential equation along the shock) is not explicitly enforced, which means that there could be errors in the numerical adjoint solution as the ones reported in [10] in 1D cases.The numerical results presented in this work and in [15] do not seem to support this concern, but a detailed grid convergence study would need to be performed to verify this assertion.Results in 1D indicate that for numerical adjoint results to converge to the analytical solution as the mesh is refined there must be consistency between the flow and adjoint calculations regarding the level of numerical smoothing (i.e., that one actually uses the same scheme with the same dissipation levels in both cases) and, more critically, that the level of dissipation increases with mesh refinement in such a way that the number of points across the shock increases while at the same time the overall width of the shock decreases.This is relatively simple to achieve n 1D, but the extension to 2D needs to be addressed with care.We hope to come back to these issues in the future.
computations reported in the paper have been carried out with the SU2 code, an open source platform developed and maintained by the SU2 Foundation.

Figure 1 .
Figure 1.Scheme of shock location and conventions.
vv py vH is the flux vector,   ˆ, , , , v ux vy p E H are the fluid's density, velocity, pressure, total energy and enthalpy, respectively,   [ a fluid perturbation U that affects both the cost function and the shock.In the perturbed flow, the new shock shape can be described in terms of a local normal deformation

Figure 2 .
Figure 2. Scheme of shock perturbation the Euler equations (2)) defined in  \ , where  denotes the fluid domain, and the shock adjoint  s , the Lagrange multiplier for the RH conditions (3), which is defined along  .With the aid of the adjoint fields, the cost function is reformulated as derivative along the shock[16],      the shock, and   b

Figure 3 .
Figure 3. Details of shock perturbation and computation of the shock displacement term

Figure 4 .
Figure 4. Close-up of the baseline NACA0012 mesh.

5 Frame 001 Figure 5 .
Figure 5. Mach contours and cutting line for flow past a NACA0012 airfoil at M = 1.5 and α = 0°.Plots of the adjoint variables are shown in Figure 6.Here, and in what follows, the adjoint variables are non-dimensional.Dimensions can be restored by multiplying the variables by the appropriate powers of   and  v , viz.

Figure 6 .
Figure 6.Flow past a NACA0012 airfoil at M = 1.5 and α = 0°.Contour lines for the adjoint x-momentum variable ψ2 (a) and y-momentum variable ψ3 (b).The bow and fishtail shocks are indicated for reference.

Figure 7 .
Figure 7. Flow past a NACA0012 airfoil at M = 1.5 and α = 0°.Contour plot for the adjoint y-momentum variable ψ3 near the bow shock.The bow shock and the cutting line are indicated for reference.

Figure 8 .
Figure 8. Flow past a NACA0012 airfoil at M = 1.5 and α = 0°.Adjoint variables along the cutting line indicated in Figure 7. Mach number and tangential velocity are also shown for reference

Figure 9 .
Figure 9. Flow past a NACA0012 airfoil at M = 1.5 and α = 0°.Adjoint normal derivatives and shock relations across the bow shock.

Figure 10 .
Figure 10.Mach contours and cutting line for flow past a NACA0012 airfoil at M = 0.8 and α = 1.25°.

Figure 11 .
Figure 11.Flow past a NACA0012 airfoil at M = 0.8 and α = 1.25°.Contour lines for the x-momentum variable ψ2 (a) and y-momentum variable ψ3 (b).The shock position is indicated for reference.

Figure 12 .
Figure 12.Flow past a NACA0012 airfoil at M = 0.8 and α = 1.25°.Contour plot for the adjoint y-momentum variable ψ3 near the shock.The shock and the cutting line are indicated for reference.

Figure 13 .Figure 14 .
Figure 13.Adjoint variables along the cutting line indicated in Figure 12.Mach number and tangential velocity are also shown for reference.

  11 MM 11 MM 11 MM
and, thus, the resulting shock is only normal at its foot.We will now consider a special incidence angle   , which should result in a normal shock to some distance above the surface.The value of   1 M is not known a priori and was obtained in[17] through an iterative process by varying the angle of attack while holding  M constant.The final value was however not disclosed, so after some numerical experiments we have found that α ≈ 5.794° results in   as illustrated in Figure15.

Figure 16 . 4 MFigure 17 .
Figure 16.Adjoint variables along the cutting line indicated in Figure 15.Mach number and tangential velocity are also shown for reference.

Acknowledgments:
The research described in this paper has been supported by INTA and the Ministry of Defence of Spain under the grants Termofluidodinámica (IGB99001) and IDATEC (IGB21001).The numerical s where (15)s generically non-zero on either side of a shock, so the only way the first equation in(15)can hold for arbitrary values of U that are also independent on either side of the shock is if