On the Characteristic Structure of the Adjoint Euler Equations and the Analytic Adjoint Solution of Supersonic Inviscid Flows

1. Introduction

This paper considers the characteristic structure of the adjoint Euler equations [1] with particular emphasis in supersonic flow, where the adjoint approach has been applied to a wide variety of problems ranging from general-purpose aerodynamic shape design [2,3], sonic boom mitigation [4,5,6,7], flow control [9], anisotropic grid adaptation [9], etc. For simplicity, the analysis is restricted to two dimensional flows. Steady inviscid compressible flow in two dimensions obeys the Euler equations

$$R(U)=\nabla \cdot \overrightarrow{F}(U)={\partial }_x{F}_x+{\partial }_y{F}_y=0$$

(1)

where

$$U=\left(\begin{array}{c}\rho \\ \rho u\\ \rho v\\ \rho E\end{array}\right),\hspace{0.33em}\hspace{0.33em}{F}_x=\left(\begin{array}{c}\rho u\\ \rho u^2+p\\ \rho uv\\ \rho uH\end{array}\right),\hspace{0.33em}\hspace{0.33em}{F}_y=\left(\begin{array}{c}\rho v\\ \rho vu\\ \rho v^2+p\\ \rho vH\end{array}\right)$$

(2)

In eq. (2), $(u,v)$ are the Cartesian components of the velocity, $\rho$ is the density, p is the pressure and E and H are the total energy and enthalpy, respectively. For a perfect gas,

$$p=(\gamma -1)\rho \left(E-{\displaystyle \frac{1}{2}}{q}^2\right)$$

(3)

and

$$\rho H=\rho E+p$$

(4)

where $q=\sqrt{u^2+v^2}$ is the modulus of the velocity and $\gamma$ is the ratio of specific heats. The flow equations (1) can be written in quasi-linear form

$$R(U)=({\partial }_U{F}_x,{\partial }_U{F}_y)\cdot ({\partial }_{x}U,{\partial }_{y}U)=0$$

(5)

where $({\partial }_U{F}_x,{\partial }_U{F}_y)$ are the inviscid flux Jacobian matrices. Eq. (5) is the starting point for the characteristic analysis of the Euler equations as we will show momentarily.

The derivation of the adjoint Euler equations is quite standard (see, for example, [1,10,11,12]) and will not be repeated here. For definiteness, we will focus on an external aerodynamics problem consisting on an airfoil with profile S immersed in a flow with incidence ${\theta }_{\infty }$. For this case, we wish to compute the sensitivities of a cost function measuring aerodynamic drag given as the integral of the pressure along a wall boundary S

$${\displaystyle \frac{1}{\kappa }}{\displaystyle {\int }_{S}p\widehat{n}\cdot \overrightarrow{d}ds}$$

(6)

where $\kappa ={\rho }_{\infty }{\left|{\overrightarrow{v}}_{\infty }\right|}^2\mathcal{l}/2$ is a normalization constant, $\mathcal{l}$ is a reference length – the chord length of the airfoil, $\widehat{n}=(n_x,n_y)$ is the wall unit normal vector and $\overrightarrow{d}=(\cos {\theta }_{\infty },\sin {\theta }_{\infty })$ lies along the inflow direction. The adjoint equations are

$$-{({\partial }_U{F}_x)}^T{\partial }_x\psi -{({\partial }_U{F}_y)}^T{\partial }_y\psi =0$$

(7)

where $\psi ={({\psi }_1,{\psi }_2,{\psi }_3,{\psi }_4)}^T$ are the adjoint variables. The adjoint boundary conditions are chosen to eliminate the integral of ${\psi }^T({\partial }_U{F}_x,{\partial }_U{F}_y)\cdot (n_x,n_y)\delta U$ on the boundaries. This procedure results in dual characteristic boundary conditions (obtained from a locally one-dimensional characteristic decomposition) in the far-field [1,10,11], while at the wall the adjoint variables obey

$$n_x{\psi }_2+n_y{\psi }_3=(n_x\cos {\theta }_{\infty }+n_y\sin {\theta }_{\infty })/\kappa$$

(8)

The characteristic structure of eq. (5) is well known and depends on the local Mach number M [13]. For subsonic flow, there is only one family of characteristic lines (the streamlines), while for supersonic flows there are two additional families of characteristics, the Mach lines that are inclined at an angle $\pm {\sin }^{-1}(1/M)$ to the local flow direction. Along characteristics, the partial differential equations (PDE) that describe the flow reduce to ordinary differential equations (ODE) called compatibility equations. In certain circumstances, notably in the case of two-dimensional irrotational supersonic flow, the compatibility conditions can be integrated and are further reduced to algebraic equations that hold only along the characteristic lines. The information carried by the compatibility conditions can be used to simplify the computation of supersonic flows, one notable practical application being the design of supersonic nozzles [13].

The adjoint equations (7) share the same characteristic structure as the flow equations but with the sign of each characteristic velocity reversed so that the characteristic information travels in the opposite direction [1,12]. As in the previous case, the adjoint variables obey compatibility conditions along the characteristics (see e.g. [14] and references therein). In the recent works [15,16], the flow and adjoint compatibility conditions were obtained and checked on several numerical adjoint solutions.

The adjoint equations are linear and, for supersonic flows, hyperbolic. As a result of this, the behavior of the analytic supersonic adjoint solutions is strongly constrained by the characteristic structure of the equations. On the one hand, and barring special cases where cost functions are defined along certain lines, not necessarily characteristic, such as in [4], the adjoint variables can only have discontinuities along characteristic lines. The jumps in the adjoint variables along these lines are constrained by jump conditions as will be shown in section 2. Additionally, as in supersonic flow information travels along characteristic lines, supersonic adjoint solutions tend to follow characteristic lines emanating from significant features of the flow or the geometry such as leading and trailing edges, shock feet, etc. as shown in, e.g., [14,15,16,17,18,1920] and in section 3. This is significant as it is precisely in those regions, where the amplitude of the adjoint solution is usually highest, that the cost function is most sensitive to perturbations.

In this paper, we review the characteristic structure of the adjoint Euler equations and obtain the corresponding compatibility and jump conditions along and across characteristic lines and show how this information can be used to obtain quantitative predictions for the adjoint variables by considering supersonic past a diamond airfoil, for which a closed-form analytic adjoint solution, valid in the near wall region, is presented and checked against a numerical adjoint solution on a very fine grid.

2. Characteristic Structure of the Adjoint Euler Equations

This section contains a brief analysis of the characteristic structure of the flow and adjoint Euler equations. We also refer the reader to the recent works [15,16], where a similar analysis is carried out.

2.1. Compatibility Conditions: Reduction of the PDE to an ODE

Both the 2D steady flow and adjoint Euler equations (5) and (7) are written in quasilinear form,

$$A{\partial }_{x}u+B{\partial }_{y}u=0$$

(9)

For a general quasilinear system (9), its characteristic structure is determined by the solution of the eigenvalue problem [21]

$$\det (B-\lambda A)=0$$

(10)

Notice that, since $\det (B-\lambda A)=\det (B^T-\lambda A^T)$, it is clear that the characteristic structure of the flow and adjoint equations is identical.

If $\lambda$ is a real solution of eq. (10), then the planar curve with local slope

$${\displaystyle \frac{dy}{dx}}=\lambda$$

(11)

is a characteristic curve of the above quasilinear system. Across characteristics, both the solution and its first derivatives may be discontinuous. Furthermore, the original system of equations can be reduced to an ODE along characteristics. Such ODES, one per characteristic curve, are called compatibility conditions, and can be obtained as follows.

Suppose $\lambda$ is a real eigenvalue of $\det (B-\lambda A)=0$. The matrix $B-\lambda A$ is singular, so there exists a left eigenvector l_λ such that

$$l_{\lambda }^T(B-\lambda A)=0$$

(12)

Multiplying eq. (9) on the left by l_λ and using the properties of the left eigenvector yields a scalar equation

$$\begin{array}{c}{l}_{\lambda }^{T}A{\partial }_{x}u+l_{\lambda }^{T}B{\partial }_{y}u=\\ l_{\lambda }^{T}A{\partial }_{x}u+\lambda l_{\lambda }^{T}A{\partial }_{y}u=\\ l_{\lambda }^{T}A({\partial }_{x}u+\lambda {\partial }_{y}u)=0\end{array}$$

(13)

In eq. (13),

$${\partial }_{x}u+\lambda {\partial }_{y}u$$

(14)

is (proportional to) the (tangent) derivative along the characteristic curve $(x,y(x))$ with local slope λ = dy/dx, so the scalar ODE (13) reduces to

$$l_{\lambda }^{T}A{\displaystyle \frac{du}{dx}}=0$$

(15)

Eq. (15) is the compatibility condition obeyed along the characteristic line by the solution to the PDE. There is one such compatibility condition along each of the families of characteristics. Besides, if the compatibility condition is integrable, it can be rearranged into the generic form

$${\displaystyle \frac{d{R}_{\lambda }}{dx}}=0$$

(16)

where $R_{\lambda }$, which is conserved along the characteristics λ = dy/dx, is called a Riemann invariant.

All the above translates immediately to the Euler equations. The adjoint equations, however, are transposed, so a bit of additional work is required. To obtain the adjoint compatibility conditions we consider again the eigenvalue problem $\det (B-\lambda A)=0$, where $(A,B)=({\partial }_U{F}_x,{\partial }_U{F}_y)$ . For real $\lambda$ one can define, in addition to the left eigenvector, a right eigenvector $r_{\lambda }$ such that

$$(B-\lambda A)r_{\lambda }=0$$

(17)

Multiplying the adjoint equation on the left by $r_{\lambda }$ and using $B{r}_{\lambda }=\lambda A{r}_{\lambda }$ yields the adjoint compatibility condition

$$r_{\lambda }^T(A^T{\partial }_x\psi +B^T{\partial }_y\psi )=r_{\lambda }^T{A}^T({\partial }_x\psi +\lambda {\partial }_y\psi )=r_{\lambda }^T{A}^T{\displaystyle \frac{d\psi }{dx}}=0$$

(18)

Under generic circumstances, it is unlikely that the adjoint compatibility conditions can be integrated. However, if the coefficients $r_{\lambda }^T{A}^T$ are constant along a certain characteristic line or in regions of constant flow, the quantities $r_{\lambda }^T{A}^T\psi$ (that we will call Riemann functions in what follows) are Riemann invariants for the adjoint equations along λ characteristics.

2.2. Discontinuities Across Characteristics: Jump Conditions

Across characteristics, both the solution and its first derivatives may be discontinuous [21], but the structure of the equations impose matching or jump conditions on the solution and its derivatives. Let us consider a characteristic line associated to the real eigenvalue$\lambda$. We assume that the Jacobian matrices (A,B) are continuous (the case of discontinuous Jacobians corresponds to shocks or slip lines that have a singular treatment of their own [22,23,24] and will not be further considered here) and distinguish two cases depending on whether $\psi$ is continuous or not.

If $\psi$ is continuous, the tangent derivatives across the characteristic line are also continuous and, by decomposing the equation in a local tangent and normal frame relative to the characteristic and taking differences across the characteristic line we get

$${\widehat{n}}_{\lambda }\cdot {(A,B)}^T{\left[{\partial }_n\psi \right]}_{\lambda }=0$$

(19)

where ${\widehat{n}}_{\lambda }\sim (-\lambda ,1)$ is the normal vector to the characteristic line.

On the other hand, if $\psi$ is discontinuous across the characteristic line, its behavior across the discontinuity is given by the following jump conditions [21]

$${\widehat{n}}_{\lambda }\cdot {(A,B)}^T{\left[\psi \right]}_{\lambda }=0$$

(20)

The proof of eq. (20) is identical to that of the Rankine-Hugoniot conditions for shocks, involving the integration of the adjoint equation against a test function on a suitably chosen domain divided in two by the line of discontinuity.

Eq. (20) is equivalent to ${(B-\lambda A)}^T{\left[\psi \right]}_{\lambda }=0$, so the jump in the adjoint variables across a characteristic line must be a left eigenvector of $B-\lambda A$.

Notice that, in either eqs. (19) and (20), the coefficient matrix of the linear system ${\widehat{n}}_{\lambda }\cdot {(A,B)}^T\sim {(-\lambda A+B)}^T$ is not invertible along the characteristic, so the above equations admit non-trivial solutions. Conversely, if the line is not characteristic, then neither eqs. (19) nor (20) admit non-trivial solutions.

2.3. Relation of Jumps Across Characteristic Lines and Compatibility Conditions

The jump conditions give a series of linear combinations of adjoint variables that are continuous across the characteristic line. Of course, not all of these are linearly independent (otherwise, the adjoint variables would be continuous across the characteristic). In fact, ${\widehat{n}}_{\lambda }\cdot {(A,B)}^T$ has rank n − k, where n is the number of equations and k is the multiplicity of the eigenvalue. This can be understood with the following result, which also sheds light onto the jump conditions. Pick a characteristic line with eigenvalue ${\lambda }_i$ and multiplicity k_i. The corresponding jump conditions are

$$({\lambda }_i{A}^T-B^T){[\psi ]}_{{\lambda }_i}=0$$

(21)

Now multiply (21) on the left by the right eigenvector corresponding to one of the n − k_i remaining eigenvalues,

$$r_{{\lambda }_j}^T({\lambda }_i{A}^T-B^T){[\psi ]}_{{\lambda }_i}=0=({\lambda }_i-{\lambda }_j)r_{{\lambda }_j}^T{A}^T{[\psi ]}_{{\lambda }_i}$$

(22)

Now, since ${\lambda }_i\ne {\lambda }_j$, the jump conditions can be written as

$$r_{{\lambda }_j}^T{A}^T{\left[\psi \right]}_{{\lambda }_i}=0$$

(23)

for j ≠ i. Now recall that along a characteristic j, the adjoint variables obey the compatibility condition $r_{{\lambda }_j}^T{A}^T\delta \psi =0$, from where we can define the Riemann functions $r_{{\lambda }_j}^T{A}^T\psi$ (which turn out to be Riemann invariants along the j characteristic under particular circumstances). We thus see that the jump conditions across one characteristic line are equivalent to the statement that the Riemann functions associated to characteristic lines that cross the given characteristic are continuous.

2.4. Characteristic Lines and Compatibility Conditions for the Adjoint Euler Equations

For the 2D Euler equations, the characteristic equation (10) has the following solutions [13]:

$${\lambda }_0={\displaystyle \frac{v}{u}}$$

(24)

(double) and

$${\lambda }_{\pm }={\displaystyle \frac{uv\pm c^2\alpha }{u^2-c^2}}={\displaystyle \frac{v\alpha \pm u}{u\alpha \mp v}}$$

(25)

where $c=\sqrt{\gamma p/\rho }$ is the speed of sound, $\alpha =\sqrt{M^2-1}$ and $M=q/c$ is the local Mach number. The characteristic lines corresponding to λ₀ are the flow streamlines. The associated compatibility conditions (15) can be arranged in the following form

$$\begin{array}{c}{K}_H=\delta H=0\\ K_S=c^{-2}\delta p-\delta \rho =0\end{array}$$

(26)

(here $\delta ={\partial }_x+\lambda {\partial }_y$ represents variation along the characteristic line) which simply say that the enthalpy and entropy are constant along streamlines (even though they might vary from streamline to streamline).

On the other hand, ${\lambda }_{\pm }$ are only real for supersonic flow, in which case the associated characteristic lines $C_{\pm }$ are, respectively, left (upper) and right (lower) -running Mach lines inclined at an angle $\pm \mu =\pm {\sin }^{-1}(1/M)$ relative to the local streamline. The corresponding compatibility conditions are

$$K_{\pm }=\delta p\mp {\displaystyle \frac{\rho }{\alpha }}(v\delta u-u\delta v)$$

(27)

For irrotational flow, eq. (27) is equivalent to

$$\delta \theta =\pm \alpha {\displaystyle \frac{\delta q}{q}}$$

(28)

where $\theta$ is the angle that the local streamline makes with the x axis. Eq. (28) can be integrated to give the Riemann invariants

$$R_{\pm }=\theta \mp \nu (M)$$

(29)

where $\nu (M)$ is the Prandtl-Meyer function [13].

For the adjoint equations, we have the exact same eigenvalues and characteristic lines, and the associated compatibility conditions $r_{\lambda }^T{A}^T\delta \psi$ can be obtained with the aid of an algebraic manipulation software (we have used Wolfgram Research’s Mathematica) and can be arranged in the following form

$$\begin{array}{l}{\lambda }_0,\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\left\{\begin{array}{c}\delta ({\psi }_1-H{\psi }_4)=0\\ \delta {\psi }_1+u\delta {\psi }_2+v\delta {\psi }_3+E_c\delta {\psi }_4=0\end{array}\right.\\ {\lambda }_{\pm },\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\delta {\psi }_1+u\delta {\psi }_2+v\delta {\psi }_3+H\delta {\psi }_4=\pm {\displaystyle \frac{1}{\alpha }}(v\delta {\psi }_2-u\delta {\psi }_3)\end{array}$$

(30)

(where $E_c=q^2/2$), which can be shown to agree with the results obtained in [15,16].

For homentropic flows, the ${\lambda }_{\pm }$ compatibility conditions can be written as

$${\lambda }_{\pm },\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\delta ({\psi }_1+u{\psi }_2+v{\psi }_3+H{\psi }_4)=\pm {\displaystyle \frac{1}{\rho \alpha }}\delta (\rho v{\psi }_2-\rho u{\psi }_3)$$

(31)

The first compatibility condition in (30) can be immediately integrated, yielding the adjoint Riemann invariant ${\psi }_1-H{\psi }_4$, which is known to be constant along streamlines (and everywhere for cost functions that only depend on the pressure [12]). Furthermore, in supersonic regions where the flow variables are constant, the above compatibility conditions give rise to 4 adjoint Riemann invariants ${r_{\lambda }}^T{A}^T\psi$

$$\begin{array}{l}{R}_1^{\psi }={\psi }_1-H{\psi }_4\\ R_2^{\psi }={\psi }_1+u{\psi }_2+v{\psi }_3+E_c{\psi }_4\end{array}$$

(32)

along streamlines,

$$R_{-}^{\psi }={\psi }_1+(u+v/\alpha ){\psi }_2+(v-u/\alpha ){\psi }_3+H{\psi }_4$$

(33)

along right-running characteristics C₋ (with eigenvalue ${\lambda }_{-}$) and

$$R_{+}^{\psi }={\psi }_1+(u-v/\alpha ){\psi }_2+(v+u/\alpha ){\psi }_3+H{\psi }_4$$

(34)

along left-running characteristics C₊ (with eigenvalue ${\lambda }_{+}$).

Even when they are not truly invariant, these Riemann functions have an additional interpretation as the adjoint variables corresponding to the flow compatibility conditions since

$$\begin{array}{l}{\psi }^T\nabla \cdot \overrightarrow{F}={\psi }^{T}A{\partial }_{x}U+{\psi }^{T}B{\partial }_{y}U={\displaystyle \frac{2u}{q^2}}{R}_2^{\psi }(\rho K_H-H{K}_S)+{\displaystyle \frac{u}{q^2}}{R}_1^{\psi }((H-E_c)K_S-\rho K_H)\\ +{\displaystyle \frac{\alpha \left(u\alpha +v\right)}{2q^2}}{R}_{-}^{\psi }{K}_{-}-{\displaystyle \frac{\alpha \left(v-u\alpha \right)}{2q^2}}{R}_{+}^{\psi }{K}_{+}\end{array}$$

(35)

Hence, $R_{1,2}^{\psi }$ are related to perturbations along streamlines, while $R_{\pm }^{\psi }$ are related to perturbations along Mach lines. We will use this insight in the next section when we analyze a particular example.

One last application of the Riemann functions concerns the jump conditions across characteristics. Across characteristic lines the adjoint solution can jump. When it does, the jumps ${\left[\psi \right]}_{\lambda }$ are subject to the conditions (20). The matrix ${\widehat{n}}_{\lambda }\cdot {(A,B)}^T$ has rank 2 for streamlines and rank 3 for Mach lines. So for a streamline, there are 2 independent jump conditions, which correspond to the continuity of the Riemann functions $R_{\pm }^{\psi }$, while for a Mach line $C_{+}$ (resp. $C_{-}$) the 3 linearly independent jump conditions can be written in terms of the Riemann functions $R_1^{\psi }$ and $R_2^{\psi }$ and $R_{-}^{\psi }$ (resp. $R_{+}^{\psi }$), which is the Riemann function corresponding to the other Mach line.

3. Application of the Jump Conditions to a Supersonic case: Flow Past a Diamond Airfoil

The method of characteristics can be used to solve for the flowfield in the case of steady, supersonic flow, and can be applied to the design of supersonic nozzles for 2D shock-free, isentropic flow [13], since in that case the non-linear flow equations reduce to algebraic equations along the characteristic lines. Since adjoint compatibility conditions cannot be integrated in general, there is little hope that a similar approach can be used with the adjoint equations. However, it turns out that the characteristic structure of the equations can be used to obtain analytic results in particular cases. One such example is given by the supersonic flow past a double wedge (diamond shaped) airfoil.

We consider inviscid supersonic flow over a symmetric diamond airfoil, with chord length of one, and a thickness of 0.06 (resulting in a half-angle of τ = 6.85 deg). The free-stream Mach number is M_∞ = 2 and the angle of attack (AOA) is 0 degrees. The flow contains two sets of wedge-shaped oblique shocks attached to the leading and trailing edges, and two expansion fans emanating at the mid-chord vertices (Figure 1).

The flow solution can be computed exactly using shock-expansion theory [13]. The flow is parallel to the free stream upstream of the leading edge shock and downstream of the trailing edge “fishtail” shock. Between the leading edge shock and the leading Mach line of the expansion fan, the flow is parallel to the front segment of the airfoil, with a Mach number M₁ < M_∞ that depends on the free stream Mach number and the wedge angle τ, which also determine the shock inclination. The flow then turns across the expansion fan (with limit Mach numbers M₁ and M₂) so that behind the trailing Mach line of the fan the flow is parallel to the rear segment of the airfoil until it reaches the fishtail shock. The Mach number along the rear section M₂ > M_∞ is determined by M₁ and the turning angle 2τ. Finally, the wedge angle, M₂ and M_∞ determine the fishtail shock inclination. From shock-expansion theory, M₁ ≈ 1.755 and M₂ ≈ 2.254, which are in very good agreement with the numerical solution shown in Figure 2, computed with DLR’s Tau solver [25] on an unstructured mesh with over 8.2 million nodes and 5300 nodes on the airfoil profile, with an outer freestream boundary domain of around 50 chord lengths from the geometry.

The corresponding drag-based adjoint solution computed with Tau’s discrete adjoint solver is shown in Figure 3. The adjoint variables are non-dimensionalized relative to Tau’s reference values $({\rho }_{\infty },p_{\infty },q_{ref}=\sqrt{p_{\infty }/{\rho }_{\infty }})$.

Figure 3 depicts the numerical adjoint solution along with the shocks and several notable characteristic lines (various streamlines, which are the nearly horizontal lines, as well as several Mach lines –the diagonal lines–, including the limits of the expansion fan and two Mach lines emanating from the trailing edge and running diagonally towards the upcoming flow). We see that the adjoint solution clearly follows the characteristic structure of the flow in a pattern that somehow mirrors that of the primal flow. The adjoint solution vanishes downstream of the two Mach lines impinging on the trailing edge since no perturbation past those lines can affect the flow about the airfoil. The solution is essentially concentrated along strips limited by Mach lines emanating from the trailing edge, the midchord vertices and the leading edge, and the adjoint solution is discontinuous along these characteristic lines. There is also a weak horizontal strip along the incoming stagnation streamline upstream of the forward shock, in agreement with the general structure described in [19]. Finally, the solution along the wall is piece-wise constant (see Figure 4), and it is actually possible to use the jump conditions (20) and the wall boundary condition to predict the values of the adjoint variables.

In this particular case, we start from the trailing edge and focus on the upper side. The right-moving Mach line emanating from the trailing edge separates two regions where the flow and the adjoint solutions are constant. Hence, the Riemann functions are piecewise constant, and the only one that can jump is the one associated to that characteristic line. The other 3 are continuous across the line and, being constant on either side, maintain the zero value that they have downstream of the Mach line. Hence, the adjoint variables upstream of the Mach line obey the equations

$$\begin{array}{l}{\psi }_1-H{\psi }_4=0\\ {\psi }_1+u{\psi }_2+v{\psi }_3+E_c{\psi }_4=0\\ {\psi }_1+(u-v/\alpha ){\psi }_2+(v+u/\alpha ){\psi }_3+H{\psi }_4=0\end{array}$$

(36)

as well as the wall boundary condition eq. (8), which on the rearmost upper segment is simply ${\psi }_2+\cot (\tau ){\psi }_3={c_{\infty }}^{-1}$ . It is now easy to see that the above equations yield for the adjoint variables in the rearmost segment of the wall the following result

$$\left(\begin{array}{c}{\psi }_1\\ {\psi }_2\\ {\psi }_3\\ {\psi }_4\end{array}\right)={\displaystyle \frac{\sin \tau }{\kappa {\alpha }_2}}\left(\begin{array}{c}-{\displaystyle \frac{\left(\gamma -1\right)M_2}{c_2}}{H}_2\\ \left(1+\left(\gamma -1\right){M_2}^2\right)\cos \tau +{\alpha }_2\sin \tau \\ {\alpha }_2\cos \tau -\left(1+\left(\gamma -1\right){M_2}^2\right)\sin \tau \\ -{\displaystyle \frac{\left(\gamma -1\right)M_2}{c_2}}\end{array}\right)$$

(37)

where ${\alpha }_2=\sqrt{{M_2}^2-1}$ and c₂ and H₂ are the speed of sound and total enthalpy at the rear part of the airfoil. Eq. (37) also holds on the lower side with the opposite sign for ${\psi }_3$. This analytic solution is compared with the numerical solution in Figure 4, showing excellent agreement.

We see in Figure 4 that the adjoint solution along the wall is piecewise constant, with the values corresponding to the rear segment of the airfoil given by eq. (37). Away from the wall, the adjoint solution exhibits a fairly simple structure, at least in the immediate proximity of the wall, as can be seen in Figure 5 where the adjoint solution along a streamline is depicted and compared with the solution along the wall. The solution is again piecewise constant, and 5 plateaus can be clearly identified, which are separated by jumps across characteristic lines. The adjoint variable is zero (V) beyond the right-running Mach line emanating from the trailing edge (7). Then it stays constant (IV) until the fan (III), delimited by Mach lines (3) and (5), where it has a continuous variation. It achieves a second plateau (II), which does not appear in the wall solution, and then jumps again across the right-running Mach line (2) emanating from the vertex of the fan [14,17]. Finally, it jumps again across the right-running Mach line emanating from the leading edge (1). Notice that the plateau values I and IV agree fairly well with the wall values, which means that the adjoint solution is roughly constant throughout the corresponding regions.

It is also interesting to note that the Riemann function $R_{+}^{\psi }$ (34), associated to left-running characteristics, is constant (and in fact vanishes) throughout the upper side of the airfoil. In fact, it is only non-zero in the region of the pressure side of the airfoil limited by the left-running characteristics emanating from the leading and trailing edges – see Figure 6. The reason for this behavior can be found in the interpretation of the adjoint Riemann functions as adjoint variables to the flow compatibility conditions. $R_{+}^{\psi }$ is the adjoint variable associated to left-running characteristics, and it is only non-zero in the region of the fluid domain where perturbations carried by left-running characteristics can reach the airfoil and, thus have an impact on the cost function.

We can use this information to extend the analytic solution to the forward part of the airfoil. We start by recalling Giles and Pierce’s Green’s function approach [22], which allows to write any 2D inviscid adjoint solution in the generic form

$$I_1\left(\begin{array}{c}1+\frac{\gamma -1}{2}{M}^2\\ -{\displaystyle \frac{u\left(1+\left(\gamma -1\right)M^2\right)}{q^2}}\\ -{\displaystyle \frac{v\left(1+\left(\gamma -1\right)M^2\right)}{q^2}}\\ {\displaystyle \frac{1+\frac{\gamma -1}{2}{M}^2}{H}}\end{array}\right)+I_2\left(\begin{array}{c}0\\ -{\displaystyle \frac{v}{q^2}}\\ {\displaystyle \frac{u}{q^2}}\\ 0\end{array}\right)+I_3\left(\begin{array}{c}-H\\ 0\\ 0\\ 1\end{array}\right)+I_4\left(\begin{array}{c}-{\displaystyle \frac{\gamma }{2}}{p}_t{M}^2\\ \gamma p_t{M}^2{\displaystyle \frac{u}{q^2}}\\ \gamma p_t{M}^2{\displaystyle \frac{v}{q^2}}\\ -{\displaystyle \frac{\gamma }{2H}}{p}_t{M}^2\end{array}\right)$$

(38)

where p_t is the total pressure (which is constant between the leading and trailing shocks) and I_j, at any given point with coordinates $\overrightarrow{x}$, are the linearized cost functions corresponding to 4 linearly independent point source perturbations (mass, force normal to the flow direction, enthalpy and total pressure) [12] inserted at $\overrightarrow{x}$. For cost functions that only depend on pressure, I₃ = 0 (total enthalpy perturbation at constant pressure), while $R_{+}^{\psi }=0$ requires that $I_2=-\sqrt{M^2-1}{I}_1$. Using this in eq. (38) and rearranging, yields

$${\displaystyle \frac{I_2}{Mc}}\left(\begin{array}{c}-{\displaystyle \frac{\left(\gamma -1\right)M}{c\alpha }}H\\ {\displaystyle \frac{1+\left(\gamma -1\right)M^2}{\alpha }}\cos \theta -\sin \theta \\ {\displaystyle \frac{1+\left(\gamma -1\right)M^2}{\alpha }}\sin \theta +\cos \theta \\ -{\displaystyle \frac{\left(\gamma -1\right)M}{c\alpha }}\end{array}\right)+I_4\left(\begin{array}{c}-{\displaystyle \frac{\gamma }{2}}{p}_t{M}^2\\ {\displaystyle \frac{\gamma p_{t}M}{c}}\cos \theta \\ {\displaystyle \frac{\gamma p_{t}M}{c}}\sin \theta \\ -{\displaystyle \frac{\gamma }{2H}}{p}_t{M}^2\end{array}\right)$$

(39)

This is the general analytic solution. It depends on two functions $I_2$ and $I_4$, whose values, for the numerical solution, are depicted in Figure 7.

In the extreme zones I and IV, $I_2$ is constant and its value is related to the adjoint wall boundary condition

$$I_2=-v{\psi }_2+u{\psi }_3=Mc({\widehat{n}}_x{\psi }_2+{\widehat{n}}_y{\psi }_3)=Mc{\widehat{n}}_x/\kappa$$

(40)

where $({\widehat{n}}_x,{\widehat{n}}_y)$ is the unit wall normal vector pointing towards the fluid, yielding

$$I_2^I=-M_1{c}_1\sin \tau /\kappa$$

(41)

and

$$I_2^{IV}=M_2{c}_2\sin \tau /\kappa$$

(42)

respectively. Within the expansion fan (zone III) we have non-constant, homentropic flow, so along right-running characteristics eq. (31) applies. Using also $R_1^{\psi }=0$, $R_{+}^{\psi }=0$ and $R_{\pm }^{\psi }=I_1\pm I_2/\alpha$, where $I_1={\psi }_1+u{\psi }_2+v{\psi }_3+H{\psi }_4$, we can obtain that

$$I_2^{III}(M)=\sigma \sqrt{{\displaystyle \frac{\alpha (M)}{\rho (M)}}}$$

(43)

where $\sigma$ is a constant. From eq. (43), we can obtain the (constant) value in zone II as

$$I_2^{II}=\sigma \sqrt{{\displaystyle \frac{\alpha (M_1)}{\rho (M_1)}}}$$

(44)

Lastly, in zone IV, $I_2^{IV}={\kappa }^{-1}{M}_2{c}_2\sin \tau$, which fixes $\sigma$ as

$$\sigma ={\displaystyle \frac{1}{\kappa }}\sqrt{{\displaystyle \frac{{\rho }_2}{{\alpha }_2}}}{M}_2{c}_2\sin \tau$$

(45)

We now turn to $I_4$. In zones I, II and IV, $I_4$ is constant, and actually $I_4^{IV}=0$ since $\Delta I_4=0$ across the right-running Mach line emanating from the trailing edge and $I_4=0$ downstream. We can obtain further information by using the jump conditions across the right-running characteristic emanating from the midchord vertex to relate the solutions in zones I and II. In terms of the Riemann functions, the jump conditions are

$$\begin{array}{l}\Delta {\psi }_1-H\Delta {\psi }_4=0\\ \Delta {\psi }_1+u\Delta {\psi }_2+v\Delta {\psi }_3+E_c\Delta {\psi }_4=0\\ \Delta {\psi }_1+(u-v/\alpha )\Delta {\psi }_2+(v+u/\alpha )\Delta {\psi }_3+H\Delta {\psi }_4=0\end{array}$$

(46)

The first and third equations in (46) are directly obeyed by the solutions on each side by construction, so we only have the middle condition

$$\Delta {\psi }_1+Mc\cos \theta \Delta {\psi }_2+Mc\sin \theta \Delta {\psi }_3+\frac{1}{2}{M}^2{c}^2\Delta {\psi }_4={\displaystyle \frac{\gamma }{\gamma -1}}{\displaystyle \frac{M^2{c}^2}{2H}}{p}_t\Delta I_4=0$$

(47)

which yields $\Delta I_4=0$ and, thus $I_4^I=I_4^{II}$, in agreement with Figure 7Across the expansion fan, $I_4$ changes between $I_4^I=I_4^{II}$ and $I_4^{IV}=0$. To obtain the value of $I_4$ across the fan we recall that $I_4$ is the effect on drag of a point perturbation to the stagnation pressure at constant static pressure [12]. Its value at a point $\overrightarrow{x}$ is given by the integration along the local streamline of secondary sources of mass and force normal to the local flow direction [26,27]

$$I_4(\overrightarrow{x})=-{\displaystyle \frac{1}{{\rho }_t}}\left({\displaystyle \underset{0}{\overset{\infty }{\int }}ds{\partial }_s{q}^{-2}{I}_1}+2{\displaystyle \underset{0}{\overset{\infty }{\int }}ds{q}^{-2}{\partial }_s\theta I_2}\right)$$

(48)

In eq. (48), ${\rho }_t$ is the total density and the integral is taken from $\overrightarrow{x}$ (s = 0) to the downstream farfield along the local streamline through $\overrightarrow{x}$. Let us consider a point inside the expansion fan. Firstly, the contribution to the integral vanishes downstream of the trailing Mach line of the fan. Secondly, the flow is isentropic, so the following relations hold along any streamline:

$$\begin{array}{c}{\partial }_s\log q={\displaystyle \frac{1}{M\left(1+\frac{\gamma -1}{2}{M}^2\right)}}{\partial }_{s}M\\ {\partial }_s\theta ={\displaystyle \frac{1}{M\left(1+\frac{\gamma -1}{2}{M}^2\right)}}{\partial }_{n}M\end{array}$$

(49)

Furthermore, inside the fan the flow variables are constant along left-running Mach lines, so we have

$$\cos \mu {\partial }_{s}M+\sin \mu {\partial }_{n}M=0\Rightarrow {\partial }_{n}M=-\cot \mu {\partial }_{s}M=-\alpha {\partial }_{s}M$$

(50)

(here μ is the Mach angle). Gathering eq. (48) – (50) and using $I_1=-I_2/\alpha$ we get

$$I_4^{III}(\eta )=-{\displaystyle \frac{2}{{\rho }_t}}{\displaystyle \underset{\eta }{\overset{{\eta }_2}{\int }}d\eta \frac{I_2}{q^{2}M\left(1+\frac{\gamma -1}{2}{M}^2\right)}({\alpha }^{-1}-\alpha ){\partial }_{\eta }M}$$

(51)

where $I_2$ is given by eq. (43) and we have changed the variable of integration to $\eta$, which is the angle that the local left-running Mach line makes with the x axis and ${\eta }_2$ is the corresponding value for the final Mach line of the fan. The final piece of information required to evaluate eq. (51) comes from the analytic solution for the expansion fan [28]

$$M(\eta )=\sqrt{1+{\displaystyle \frac{\gamma +1}{\gamma -1}}{\tan }^{2}z(\eta )}$$

(52)

where

$$z(\eta )=\sqrt{{\displaystyle \frac{\gamma -1}{\gamma +1}}}(\nu (M_1)+\tau -\eta +\pi /2)$$

(53)

Differentiating eq. (52) with respect to $\eta$ yields

$${\partial }_{\eta }M=-\sqrt{{\displaystyle \frac{\gamma +1}{\gamma -1}}}{\displaystyle \frac{\tan z}{M{\cos }^{2}z}}$$

(54)

Gathering all the information, we can write the analytic adjoint solutions as

$${\Psi }_I=-{\displaystyle \frac{\sin \tau }{\kappa }}\left(\begin{array}{c}-{\displaystyle \frac{\left(\gamma -1\right)M_1}{c_1{\alpha }_1}}{H}_1\\ {\displaystyle \frac{1+\left(\gamma -1\right){M_1}^2}{{\alpha }_1}}\cos \tau -\sin \tau \\ {\displaystyle \frac{1+\left(\gamma -1\right){M_1}^2}{{\alpha }_1}}\sin \tau +\cos \tau \\ -{\displaystyle \frac{\left(\gamma -1\right)M_1}{c_1{\alpha }_1}}\end{array}\right)+\gamma p_t{M}_1{I}_4^{III}({\eta }_1)\left(\begin{array}{c}-{\displaystyle \frac{M_1}{2}}\\ {\displaystyle \frac{\cos \tau }{c_1}}\\ {\displaystyle \frac{\sin \tau }{c_1}}\\ -{\displaystyle \frac{M_1}{2H_1}}\end{array}\right)$$

(55)

in zone I, where ${\eta }_1$ gives the inclination of the fan’s leading Mach line relative to the x axis,

$${\Psi }_{II}={\displaystyle \frac{\sin \tau }{\kappa }}\sqrt{{\displaystyle \frac{{\rho }_2{\alpha }_1}{{\rho }_1{\alpha }_2}}}{\displaystyle \frac{M_2{c}_2}{M_1{c}_1}}\left(\begin{array}{c}-{\displaystyle \frac{\left(\gamma -1\right)M_1}{c_1{\alpha }_1}}{H}_1\\ {\displaystyle \frac{1+\left(\gamma -1\right){M_1}^2}{{\alpha }_1}}\cos \tau -\sin \tau \\ {\displaystyle \frac{1+\left(\gamma -1\right){M_1}^2}{{\alpha }_1}}\sin \tau +\cos \tau \\ -{\displaystyle \frac{\left(\gamma -1\right)M_1}{c_1{\alpha }_1}}\end{array}\right)+\gamma p_t{M}_1{I}_4^{III}({\eta }_1)\left(\begin{array}{c}-{\displaystyle \frac{M_1}{2}}\\ {\displaystyle \frac{\cos \tau }{c_1}}\\ {\displaystyle \frac{\sin \tau }{c_1}}\\ -{\displaystyle \frac{M_1}{2H_1}}\end{array}\right)$$

(56)

in zone II,

$${\Psi }_{III}={\displaystyle \frac{\sin \tau }{\kappa }}\sqrt{{\displaystyle \frac{{\rho }_2}{{\alpha }_2}}{\displaystyle \frac{\alpha }{\rho }}}{\displaystyle \frac{M_2{c}_2}{Mc}}\left(\begin{array}{c}-{\displaystyle \frac{\left(\gamma -1\right)M}{c\alpha }}H\\ {\displaystyle \frac{1+\left(\gamma -1\right)M^2}{\alpha }}\cos \theta -\sin \theta \\ {\displaystyle \frac{1+\left(\gamma -1\right)M^2}{\alpha }}\sin \theta +\cos \theta \\ -{\displaystyle \frac{\left(\gamma -1\right)M}{c\alpha }}\end{array}\right)+I_4^{III}(\eta )\left(\begin{array}{c}-{\displaystyle \frac{\gamma }{2}}{p}_t{M}^2\\ {\displaystyle \frac{\gamma p_{t}M}{c}}\cos \theta \\ {\displaystyle \frac{\gamma p_{t}M}{c}}\sin \theta \\ -{\displaystyle \frac{\gamma }{2H}}{p}_t{M}^2\end{array}\right)$$

(57)

across the fan and

$${\Psi }_{IV}={\displaystyle \frac{\sin \tau }{\kappa }}\left(\begin{array}{c}-{\displaystyle \frac{\left(\gamma -1\right)M_2}{c_2{\alpha }_2}}{H}_2\\ {\displaystyle \frac{1+\left(\gamma -1\right){M_2}^2}{{\alpha }_2}}\cos \tau +\sin \tau \\ -{\displaystyle \frac{1+\left(\gamma -1\right){M_2}^2}{{\alpha }_2}}\sin \tau +\cos \tau \\ -{\displaystyle \frac{\left(\gamma -1\right)M_2}{c_2{\alpha }_2}}\end{array}\right)$$

(58)

in zone IV, in agreement with eq. (37). The analytic solution eq. (55) – (58) is in remarkable agreement with the numerical solution as shown in Figure 4 and Figure 8.

Solution Across the Expansion Fan

Eq. (57) gives the spatially varying solution throughout the fan. Notice that, as is clear from eq. (57), the adjoint solution, like the flow solution, remains constant along each Mach line of the fan, as can be confirmed with the numerical solution in Figure 9, which depicts I₂ and I₄ and the adjoint variables along 3 Mach lines of the fan (roughly corresponding to the leading and trailing lines of the fan and a third interior line). The corresponding analytic values for I₂ and I₄ computed with eq. (43) and (51) are also shown for comparison, showing and excellent agreement.

We see that towards the center of the fan (located at x = 0.5 in the plot), I₂ and I₄ and the adjoint variables attain a constant value along each line of the fan, the value changing from line to line. I₂ and I₄ behave differently in this regard, since the former jumps abruptly across the right-running Mach line (7), while the latter shows a smooth variation and reaches a plateau at the intersection with the limiting streamline (6). The reason for this different behavior can be explained as follows. Both I₂ and I₄ can be written in terms of the adjoint Riemann functions as

$$\begin{array}{l}{I}_2={\displaystyle \frac{\alpha }{2}}(R_{+}^{\psi }-R_{-}^{\psi })\\ I_4={\displaystyle \frac{1}{p_0}}\left({\displaystyle \frac{\gamma -1}{\gamma }}+{\displaystyle \frac{2}{\gamma M^2}}\right)R_2^{\psi }-{\displaystyle \frac{1}{\gamma p_0{M}^2}}{R}_1^{\psi }\end{array}$$

(59)

Let us now focus on the fan of the upper side of the airfoil (Figure 10). Since both $R_{+}^{\psi }$ and $R_1^{\psi }$ vanish everywhere on the upper side of the airfoil, it turns out that in the fan $I_2\sim R_{-}^{\psi }$ and $I_4\sim R_2^{\psi }$. Now, within the fan $R_{-}^{\psi }$ and $R_2^{\psi }$ are only different from zero in the shaded triangle bounded by Mach lines 3, 5 and 7. While $R_2^{\psi }$ (and, thus, I₄) is continuous across Mach line 7, $R_{-}^{\psi }$ (and, thus, I₂) is not, and its jump depends on the local value of the flow, which explains the difference between the plateau values along each line of the fan.

4. Conclusions

The behavior of the steady solutions to the adjoint Euler equations is severely constrained by the characteristic structure of the equations, particularly in supersonic flows, where the equations are hyperbolic and the solutions display distinctive traits along certain significant characteristic lines.

In this work, we have analyzed the characteristic structure of the adjoint Euler equations in two dimensions. The eigenvalues and characteristic lines are the same as for the Euler equations, and compatibility conditions can be derived that constrain the evolution of the adjoint variables along the characteristics. At each point, these compatibility conditions are equivalent to the adjoint Euler equations. Additionally, adjoint solutions can be discontinuous along characteristic lines, and the jumps of the adjoint variables along these lines are constrained by jump conditions that can be used to build analytic adjoint solutions in simple cases, such as the supersonic flow past a diamond airfoil considered in this work.

The interest of the results described in this work is twofold. On the one hand, they provide benchmark solutions and constrains that can be used for verification of numerical adjoint solvers. On the other hand, these insights may significantly impact methodologies used in the aerodynamic design processes and related applications in supersonic regimes, particularly as commercial supersonic flight becomes more relevant.