# Adjoint and Direct Characteristic Equations for Two-Dimensional Compressible Euler Flows

Kevin Ancourt, Jacques Peter, Olivier Atinault

► To cite this version:

Kevin Ancourt, Jacques Peter, Olivier Atinault. Adjoint and Direct Characteristic Equations for Two-Dimensional Compressible Euler Flows. Aerospace, 2023, 10 (9), pp.797. 10.3390/aerospace10090797. [https://doi.org/10.3390/aerospace10090796](https://doi.org/10.3390/aerospace10090796)

HAL Id: hal-04209946

https://hal.science/hal-04209946v1

Submitted on 18 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Article

# Adjoint and Direct Characteristic Equations for Two-Dimensional Compressible Euler Flows

Kevin Ancourt<fnref n="1" />, Jacques Peter <fnref n="1" />,* and Olivier Atinault<fnref n="2" />

**Abstract:** The method of characteristics is a classical method for gaining understanding in the solution of a partial differential equation. It has recently been applied to the adjoint equations of the 2D steady-state Euler equations and the first goal of this paper is to present a linear algebra analysis that greatly simplifies the discussion of the number of independent characteristic equations satisfied along a family of characteristic curves. This method may be applied for both the direct and the adjoint problem. Our second goal is to directly derive in conservative variables the characteristic equations of 2D compressible inviscid flows. Finally, the theoretical results are assessed for a nozzle flow with a classical scheme and its dual consistent discrete adjoint.

**Keywords:** continuous adjoint; inviscid flow; compressible flow; method of characteristics; characteristic equation; characteristic curve

## 1. Introduction

The method of characteristics is a well-known method for studying partial differential equations (PDEs). It aims to exhibit specific hypersurfaces in the input domain where the solution of the PDE of interest satisfies an ordinary differential equation (ODE). When applied to 2D steady-state inviscid compressible flows, it is known to provide a full resolution of a supersonic area based only on the knowledge of the inflow (whereas it provides partial information for a subsonic flow) [1–4]. Both a theoretical understanding and practical calculations of a variety of flows (in nozzles, along steps, along curved walls) are enabled by this technique.

In addition, discrete and continuous adjoint are now well-established methods for shape optimization [5–9] and goal-oriented mesh adaptation [10–12]. Because of these important applications, regular efforts have been devoted to the fast and safe writing of adjoint modules [13–16] and to the efficient solving of the adjoint equations [17,18]. Adjoint methods are also useful for flow control [19,20], meta-modeling [21,22], receptivity–sensitivity–stability analyses [23,24], and data assimilation [25,26].

Moving from direct to adjoint Euler equations, the classical flux Jacobians are replaced in the PDE by the opposite of their transpose. This allows the simple adaptation in 2D and 3D of the characteristic equations in the specific sense of the quasi-one-dimensional scalar propagation equations satisfied in the direction normal to a boundary [27]. From the earliest works on the continuous adjoint method, it has been proved that for the adjoint problem (a) information travels along characteristics in the opposite direction to the flow, and (b) the right eigenvectors of the local Jacobian replace the left eigenvectors in the equation definition [5,28,29]. In addition, in a recent paper, the characteristic equations (CEs), in the sense of nonlinear ODEs satisfied along curves to be defined, have been derived for the adjoint 2D steady-state Euler equations [30].

<figure><img src="image_2.png" /></figure>

### check for updates

**Citation:** Ancourt, K.; Peter, J.; Atinault, O. Adjoint and Direct Characteristic Equations for Two-Dimensional Compressible Euler Flows. *Aerospace* **2023**, *10*, 797. https://doi.org/10.3390/aerospace10090797

Academic Editor: Carlos Lozano

Received: 20 July 2023

Revised: 8 September 2023

Accepted: 10 September 2023

Published: 12 September 2023

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">¹ Département Aérodynamique, Aéroélasticité et Acoustique (DAAA), ONERA, Université Paris Saclay, 29 Avenue de la Division Leclerc, Châtillon, 92322 Paris, France; kevin.ancourt@onera.fr</footnote>

<footnote n="2">² Département Aérodynamique, Aéroélasticité et Acoustique (DAAA), ONERA, Université Paris Saclay, 8 Rue des Vertugadins, 92190 Meudon, France; olivier.atinault@onera.fr</footnote>

<footnote>* Correspondence: jacques.peter@onera.fr</footnote>
The classical direct characteristic equations (DCEs) for 2D inviscid compressible flows and the recently presented corresponding adjoint characteristic equations (ACEs) appear to be linked; in particular, the characteristic curves are the same for both problems [30]. The analogies between ACEs and DCEs are further studied in Section 2, where several properties of a generic system of equations embedding ACEs and DCEs are demonstrated. In particular, this formal linear algebra approach greatly simplifies the determination of the number of independent CEs satisfied along two families of characteristic curves.

Whereas the ACEs need to be derived from the complete 8 × 8 linear system for the Cauchy problem, the DCEs have been derived with various mechanical assumptions and also various sets of variables [3,4,31,32]. For the sake of a full understanding, the direct characteristic problem is solved in conservative variables from the complete Cauchy problem and the resulting differential forms are discussed in Section 3.

Finally, an inviscid supersonic nozzle flow is considered in Section 4, where both the DCEs and the ACEs are numerically assessed. We conclude by discussing the interest in the ACEs for code verification [33] of the continuous and discrete adjoint modules.

## 2. Common Properties of Linear Systems for the Cauchy Direct and Adjoint Problems

### 2.1. Notation, Steady-State Euler Equations

We denote by $W = (\rho, \rho u, \rho v, \rho E)$ the conservative variables, with $\rho$ being the density, $(u, v)$ the components of the velocity $\bar{U}$, $e$ the internal energy, and $E$ the total energy. A thermally and calorically perfect gas law is considered. The static pressure $p$, total enthalpy $H$, and entropy $S$ are

$$p = (\gamma - 1)\rho e = (\gamma - 1)(\rho E - 0.5\rho ||\bar{U}||^2),$$ $$H = E + \frac{p}{\rho}, \quad S = c_v \ln\left(\frac{p}{\rho\gamma}\right),$$

with a constant $\gamma = 1.4$ and a constant heat capacity at constant volume $c_v$. A subscript $i$ is used to denote the classical stagnation quantities $(\rho_i, p_i, T_i)$. The steady-state Euler equations read

$$\frac{\partial \mathbf{F}}{\partial x} + \frac{\partial \mathbf{G}}{\partial y} = 0 \qquad (1)$$

with $\mathbf{F}$ and $\mathbf{G}$ the fluxes along $x$ and $y$,

$$\mathbf{F} = \begin{pmatrix} \rho u \\ p + \rho u^2 \\ \rho u v \\ \rho u H \end{pmatrix}, \quad \mathbf{G} = \begin{pmatrix} \rho v \\ \rho u v \\ p + \rho v^2 \\ \rho v H \end{pmatrix}, \qquad (2)$$

or equivalently,

$$A(W) \frac{\partial W}{\partial x} + B(W) \frac{\partial W}{\partial y} = 0. \qquad (3)$$

The Euler flux Jacobian matrices in the $x$- and $y$-directions, $A$ and $B$, are equal to

$$A = \begin{bmatrix} 0 & 1 & 0 & 0 \\ \gamma_1 E_c - u^2 & (3 - \gamma)u & -\gamma_1 v & \gamma_1 \\ -uv & v & u & 0 \\ (\gamma_1 E_c - H)u & H - \gamma_1 u^2 & -\gamma_1 u v & \gamma u \end{bmatrix}, \quad B = \begin{bmatrix} 0 & 0 & 1 & 0 \\ -uv & v & u & 0 \\ \gamma_1 E_c - v^2 & -\gamma_1 u & (3 - \gamma)v & \gamma_1 \\ (\gamma_1 E_c - H)v & -\gamma_1 u v & H - \gamma_1 v^2 & \gamma v \end{bmatrix}$$

with $\gamma_1 = \gamma - 1$. The columns of the Jacobians are denoted by $A_1$ to $A_4$ and $B_1$ to $B_4$ such that $A = [A_1|A_2|A_3|A_4]$, $B = [B_1|B_2|B_3|B_4]$. Finally, the discrete adjoint vector of the functional output of interest is denoted as $\psi = (\psi_1, \psi_2, \psi_3, \psi_4)$.
## 2.2. Generic Cauchy Problem Embedding ACE and DCE Base Equations

Given $a$ and $b$, two neighboring points in the fluid domain, we denote by $(dx, dy) = ab$ and then $t = dy/dx$ and $\kappa = ut - v$. The increment in the conservative variables between the two points is denoted as $(d\rho, d\rho u, d\rho v, d\rho E) = (\rho^b - \rho^a, \rho u^b - \rho u^a, \rho v^b - \rho v^a, \rho E^b - \rho E^a)$ and the increment in the adjoint vector components as $(d\psi_1, d\psi_2, d\psi_3, d\psi_4) = (\psi_1^b - \psi_1^a, \psi_2^b - \psi_2^a, \psi_3^b - \psi_3^a, \psi_4^b - \psi_4^a)$. The Cauchy problem for the conservative variables aims at calculating their derivatives at point $a$ with first order in space, and the singularities of this problem are at the core of the search of the CEs. The corresponding 8 × 8 linear system associates first-order Taylor expansions and the 2D steady-state Euler equations. The complete system is

$$ \begin{bmatrix} dx & 0 & 0 & 0 & dy & 0 & 0 & 0 \\ 0 & dx & 0 & 0 & dy & 0 & 0 & 0 \\ 0 & 0 & dx & 0 & 0 & dy & 0 & 0 \\ 0 & 0 & 0 & dx & 0 & 0 & dy & 0 \\ & & & & & & & \\ A & & & & & & & B \\ & & & & & & & \end{bmatrix} \begin{bmatrix} (\partial\rho/\partial x) \\ (\partial\rho u/\partial x) \\ (\partial\rho v/\partial x) \\ (\partial\rho E/\partial x) \\ (\partial\rho/\partial y) \\ (\partial\rho u/\partial y) \\ (\partial\rho v/\partial y) \\ (\partial\rho E/\partial y) \end{bmatrix} = \begin{bmatrix} d\rho \\ d\rho u \\ d\rho v \\ d\rho E \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}. \quad (4) $$

The corresponding linear system for the derivatives of the adjoint vector reads

$$ \begin{bmatrix} dx & 0 & 0 & 0 & dy & 0 & 0 & 0 \\ 0 & dx & 0 & 0 & dy & 0 & 0 & 0 \\ 0 & 0 & dx & 0 & 0 & dy & 0 & 0 \\ 0 & 0 & 0 & dx & 0 & 0 & dy & 0 \\ & & & & & & & \\ - & A^T & & & & & & B^T \\ & & & & & & & \end{bmatrix} \begin{bmatrix} (\partial\psi_1/\partial x) \\ (\partial\psi_2/\partial x) \\ (\partial\psi_3/\partial x) \\ (\partial\psi_4/\partial x) \\ (\partial\psi_1/\partial y) \\ (\partial\psi_2/\partial y) \\ (\partial\psi_3/\partial y) \\ (\partial\psi_4/\partial y) \end{bmatrix} = \begin{bmatrix} d\psi_1 \\ d\psi_2 \\ d\psi_3 \\ d\psi_4 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}, \quad (5) $$

where $-A^T$ and $-B^T$ are the opposite and transposed Jacobians of the Euler fluxes. The straightforward counterpart for the derivatives of the discrete flow field $W$, appears in Equation (4). Several results regarding the minors of the matrices in (4) and (5), and also the number of independent CEs when $(dx, dy)$ makes the 8 × 8 matrices singular, may be demonstrated by considering a generic matrix $K$ and a generic linear system

$$ K \begin{bmatrix} (\partial \zeta_1 / \partial x) \\ (\partial \zeta_2 / \partial x) \\ (\partial \zeta_3 / \partial x) \\ (\partial \zeta_4 / \partial x) \\ (\partial \zeta_1 / \partial y) \\ (\partial \zeta_2 / \partial y) \\ (\partial \zeta_3 / \partial y) \\ (\partial \zeta_4 / \partial y) \end{bmatrix} = \begin{bmatrix} dx & 0 & 0 & 0 & dy & 0 & 0 & 0 \\ 0 & dx & 0 & 0 & dy & 0 & 0 & 0 \\ 0 & 0 & dx & 0 & 0 & dy & 0 & 0 \\ 0 & 0 & 0 & dx & 0 & 0 & dy & 0 \\ & & & & & & & \\ A & & & & & & & B \\ & & & & & & & \end{bmatrix} \begin{bmatrix} (\partial \zeta_1 / \partial x) \\ (\partial \zeta_2 / \partial x) \\ (\partial \zeta_3 / \partial x) \\ (\partial \zeta_4 / \partial x) \\ (\partial \zeta_1 / \partial y) \\ (\partial \zeta_2 / \partial y) \\ (\partial \zeta_3 / \partial y) \\ (\partial \zeta_4 / \partial y) \end{bmatrix} = \begin{bmatrix} d\zeta_1 \\ d\zeta_2 \\ d\zeta_3 \\ d\zeta_4 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}, $$

where $(A, B)$ stands either for $(-A^T, -B^T)$ or for $(A, B)$. The corresponding generic notation for $\psi$ or $W$ is $\zeta$. The columns of $A$ and $B$ are denoted as $\{A_1, A_2, A_3, A_4\}$ and $\{B_1, B_2, B_3, B_4\}$, respectively. We stress that the bold characters are, therefore, associated with the generic problem and not with a higher tensor order.
## 2.3. Coefficients of ACE and DCE as 4 × 4 Determinants

It has already been observed that the direct and adjoint problems share the same determinant, computed using linear combinations of columns, as

$$|\mathbf{K}| = |-dx\mathbf{B} + dy\mathbf{A}| = (-v \, dx + u \, dy)^2 (-v \, dx + u \, dy + c \, dl)(-v \, dx + u \, dy - c \, dl),$$

with

$$c = \sqrt{\frac{\gamma p}{\rho}}, \quad dl = \sqrt{dx^2 + dy^2},$$

using the known eigenvalues of the Euler flux Jacobian in an arbitrary direction [30]. If $|\mathbf{K}| \neq 0$, the Cauchy problem is well posed, and the solution is expressed using the minors of $\mathbf{K}$. These minors are denoted as $\mathbf{K}_{jx}^i$ and $\mathbf{K}_{jy}^i$ with indices referring to the considered variable and the direction of differentiation. For example, $(\partial\zeta_1/\partial x)$ is

$$\left(\frac{\partial \zeta_1}{\partial x}\right) = \frac{\begin{vmatrix} d\zeta_1 & 0 & 0 & 0 & dy & 0 & 0 & 0 \\ d\zeta_2 & dx & 0 & 0 & 0 & dy & 0 & 0 \\ d\zeta_3 & 0 & dx & 0 & 0 & 0 & dy & 0 \\ d\zeta_4 & 0 & 0 & dx & 0 & 0 & 0 & dy \\ | & | & | & | & | & | & | & | \\ 0 & \mathbf{A}_2 & \mathbf{A}_3 & \mathbf{A}_4 & \mathbf{B}_1 & \mathbf{B}_2 & \mathbf{B}_3 & \mathbf{B}_4 \\ | & | & | & | & | & | & | & | \end{vmatrix}}{|\mathbf{K}|} \quad (6)$$

$$\left(\frac{\partial \zeta_1}{\partial x}\right) = \frac{\mathbf{K}_{1x}^1 d\zeta_1 - \mathbf{K}_{1x}^2 d\zeta_2 + \mathbf{K}_{1x}^3 d\zeta_3 - \mathbf{K}_{1x}^4 d\zeta_4}{|\mathbf{K}|} \quad (7)$$

Conversely, the generic Cauchy problem is ill posed if and only if $|\mathbf{K}| = 0$, that is, along the $S$, $C^+$, and $C^-$ curves:

$$-v \, dx + u \, dy = 0 \qquad S \text{ streamtraces} \quad \text{(all Mach numbers)} \quad (8)$$

$$-v \, dx + u \, dy + c \, dl = 0 \qquad C^{-} \text{ characteristics} \quad \text{(supersonic flow only)} \quad (9)$$

$$-v \, dx + u \, dy - c \, dl = 0 \qquad C^{+} \text{ characteristics} \quad \text{(supersonic flow only)} \quad (10)$$

In such cases, the boundedness of $(\partial\zeta_j/\partial x)$ $(\partial\zeta_j/\partial y)$ along the characteristic curves and their expression outside these curves yields the CE. Typically, Equation (7), that is valid outside the characteristic curves, yields the following CE:

$$\mathbf{K}_{1x}^1 d\zeta_1 - \mathbf{K}_{1x}^2 d\zeta_2 + \mathbf{K}_{1x}^3 d\zeta_3 - \mathbf{K}_{1x}^4 d\zeta_4 = 0$$

along the characteristic curves. Of course, seven other CEs are derived from the boundedness of $(\partial\zeta_2/\partial x)$... $(\partial\zeta_4/\partial y)$ and we need to determine a minimal set of independent CEs for each type of characteristic curve.

To that end, the $\mathbf{K}_{jx}^i$ and $\mathbf{K}_{jy}^i$ minors are expressed as determinants of 4 × 4 matrices. From Equations (6) and (7),

$$\mathbf{K}_{1x}^1 = \begin{vmatrix} dx & 0 & 0 & 0 & dy & 0 & 0 \\ 0 & dx & 0 & 0 & 0 & dy & 0 \\ 0 & 0 & dx & 0 & 0 & 0 & dy \\ | & | & | & | & | & | & | \\ \mathbf{A}_2 & \mathbf{A}_3 & \mathbf{A}_4 & \mathbf{B}_1 & \mathbf{B}_2 & \mathbf{B}_3 & \mathbf{B}_4 \\ | & | & | & | & | & | & | \end{vmatrix}$$

is the basic formula for $\mathbf{K}_{1x}^1$. In all 32 corresponding expressions, the $dy$ terms of the first three lines are eliminated by linear combinations of columns and the determinants are then expanded along these lines. We assume that $dx \neq 0$ using $t = dy/dx$ and will further
consider this hypothesis in the following. For the sake of brevity, we shall not discuss the simple specific case where $dx = 0$.

The minors arising in the CE derived from the existence of $(\partial\zeta_1/\partial x)$ then read

$$ 
\begin{aligned}
\mathbf{K}_{1x}^1 &= dx^3 & |\mathbf{B}_1 (\mathbf{B}_2 - t\mathbf{A}_2) & (\mathbf{B}_3 - t\mathbf{A}_3) (\mathbf{B}_4 - t\mathbf{A}_4)| \\
\mathbf{K}_{1x}^2 &= -dx^2 dy & |\mathbf{A}_2 \mathbf{B}_2 & (\mathbf{B}_3 - t\mathbf{A}_3) (\mathbf{B}_4 - t\mathbf{A}_4)| \\
\mathbf{K}_{1x}^3 &= -dx^2 dy & |\mathbf{A}_3 (\mathbf{B}_2 - t\mathbf{A}_2) & \mathbf{B}_3 (\mathbf{B}_4 - t\mathbf{A}_4)| \\
\mathbf{K}_{1x}^4 &= -dx^2 dy & |\mathbf{A}_4 (\mathbf{B}_2 - t\mathbf{A}_2) & (\mathbf{B}_3 - t\mathbf{A}_3) \mathbf{B}_4|
\end{aligned}
$$

and the corresponding expressions for $(\partial\zeta_1/\partial y)$ are

$$ 
\begin{aligned}
\mathbf{K}_{1y}^1 &= -dx^3 |\mathbf{A}_1 & (\mathbf{B}_2 - t\mathbf{A}_2) & (\mathbf{B}_3 - t\mathbf{A}_3) & (\mathbf{B}_4 - t\mathbf{A}_4)| \\
\mathbf{K}_{1y}^2 &= dx^3 |\mathbf{A}_2 & \mathbf{B}_2 & (\mathbf{B}_3 - t\mathbf{A}_3) & (\mathbf{B}_4 - t\mathbf{A}_4)| \\
\mathbf{K}_{1y}^3 &= dx^3 |\mathbf{A}_3 & (\mathbf{B}_2 - t\mathbf{A}_2) & \mathbf{B}_3 & (\mathbf{B}_4 - t\mathbf{A}_4)| \\
\mathbf{K}_{1y}^4 &= dx^3 |\mathbf{A}_4 & (\mathbf{B}_2 - t\mathbf{A}_2) & (\mathbf{B}_3 - t\mathbf{A}_3) & \mathbf{B}_4|
\end{aligned}
$$

It is clear that

$$ \mathbf{K}_{1x}^2 = -t\mathbf{K}_{1y}^2 \qquad \mathbf{K}_{1x}^3 = -t\mathbf{K}_{1y}^3 \qquad \mathbf{K}_{1x}^4 = -t\mathbf{K}_{1y}^4. \qquad (11) $$

In addition, expanding $|\mathbf{K}|$ along the first line of the matrix results in

$$ |\mathbf{K}| = dx\mathbf{K}_{1x}^1 + dy\mathbf{K}_{1y}^1, $$

so that along a characteristic curve $\mathbf{K}_{1x}^1 = -t\mathbf{K}_{1y}^1$. Corresponding equations are found for the other variables $\zeta_j$ $j \ge 2$: $\mathbf{K}_{jx}^i = -t\mathbf{K}_{jy}^i$ if $i \neq j$ whatever the value of $|\mathbf{K}|$ and $\mathbf{K}_{jx}^i = -t\mathbf{K}_{jy}^i$ if $|\mathbf{K}| = 0$. If the minors are not all null, this proves that the four CEs derived from the boundedness of $(\partial\zeta_k/\partial x)$ and their counterparts derived from the boundedness of $(\partial\zeta_k/\partial y)$ are proportional, so that only one set has to be studied. At this point, the case of the S curves and the one of the $C^+$ and $C^-$ curves shall be discussed separately.

## 2.4. Number of independent CEs along the $C^+$ and $C^-$

The abstract linear algebra point of view allows us to calculate the number of independent CEs in the specific case of the $C^+$ and $C^-$ curves. For $dx \neq 0$, $\mathbf{K}$ is easily found to be equivalent to

$$ \left[ \begin{array}{cc} dx\mathbf{I} & 0 \\ \mathbf{A} & \mathbf{B} - t\mathbf{A} \end{array} \right] \quad \text{and then to} \quad \left[ \begin{array}{cc} dx\mathbf{I} & 0 \\ 0 & \mathbf{B} - t\mathbf{A} \end{array} \right]. \qquad (12) $$

When the value $t$ is the one corresponding to a $C^+$ or a $C^-$, for both the adjoint and the direct problem, the rank of $\mathbf{B} - t\mathbf{A}$ is known to be three from the eigenanalysis of the Euler equations' flux Jacobians. Using classical theorems for the ranks of block diagonal matrices and the ranks of equivalent matrices, the rank of $\mathbf{K}$ is easily proved to be seven, one less than its size. In this case, the adjugate matrix is known to have rank one. This means that all CEs (whose coefficients appear in the rows of the adjugate matrix) are proportional and reduced to only one independent equation. This result is very valuable as
the corresponding explicit calculations from the algebraic expressions of the minors are
very tedious. It has been used in Section 3.4.

2.5. Nullity of All Minors along the S Curves

The 4 × 4 determinant involved in the expression of $K_{1x}^1$ is

$$
|\mathbf{B}_1 \ (\mathbf{B}_2 - t\mathbf{A}_2) \ (\mathbf{B}_3 - t\mathbf{A}_3) \ (\mathbf{B}_4 - t\mathbf{A}_4)|. \qquad (13)
$$

Along a streamtrace $t = dy/dx = v/u$ and $(\mathbf{B} - t\mathbf{A})$ is a matrix of rank two. As all its minors of rank three are null, the determinant (13), calculated by expansion along the first column, appears to be zero. To obtain the same property for $\mathbf{K}_{1x}^2$, $\mathbf{K}_{1x}^3$, and $\mathbf{K}_{1x}^4$, a straightforward rewriting is needed. For example, the determinant in $\mathbf{K}_{1x}^2$ reads

$$
| \mathbf{A}_2 \ \mathbf{B}_2 \ (\mathbf{B}_3 - t\mathbf{A}_3) \ (\mathbf{B}_4 - t\mathbf{A}_4) | \quad \text{or possibly} \quad | \mathbf{A}_2 \ (\mathbf{B}_2 - t\mathbf{A}_2) \ (\mathbf{B}_3 - t\mathbf{A}_3) \ (\mathbf{B}_4 - t\mathbf{A}_4) |
$$

Under this latter form, it is clearly null when the rank of $(\mathbf{B} - t\mathbf{A})$ is two. Using the exact same arguments, it is possible to prove that all $K_{jx}^i$ and $K_{ix}^i$ are null when $t = v/u$.

In addition, we may change our point of view and consider (13) as a function of $t$
for fixed $W$. This function is obviously a polynomial of maximum degree three and $v/u$
is one of its roots. We, therefore, expect $(t - v/u)$ to be a factor of all $K_{jx}^i$, $K_{jy}^i$ algebraic
expressions. The explicit expressions of the coefficients for the ACEs (see [30]) and DCEs
(see Appendix A) confirm this property.

Along the S curves, **K** is found to have rank six thanks to the technique used in the previous subsection. As already stated, all $K_{jx}^i$ and $K_{jy}^i$ are zero and the adjugate matrix of **K** is the null matrix. The derivation of the explicit CE along the trajectories then involves the multiplicity two of $(-vdx + udy)$ and simplified $\bar{K}_{jx}^i$, $\bar{K}_{jy}^i$ coefficients derived by removing $(-vdx + udy)$ in the expressions of the corresponding $K_{jx}^i$, $K_{jy}^i$. The proportionality of the explicit CEs derived from the boundedness of $(\partial\zeta_k/\partial x)$ and $(\partial\zeta_k/\partial y)$ is verified but is more easily presented based on the actual formulas of the CEs—see Section 3 for the DCEs and in Section 2 of [30] for the ACEs.

3. Derivation of Characteristic Equations for 2D Inviscid Compressible Flows Using
Conservative Variables

3.1. Flow Properties, Sets of Variables, Resulting Equations for the Usual Derivation of the Characteristic Equations

The classical derivations of the DCEs use a set of Taylor expansions and mechanical equations for the variations of a set of primitive variables. These resulting equations are linear in the unknown derivatives of the primitive variables with nonlinear functions of the state variables as coefficients. The variations are generally expressed in an orthonormal frame of reference ($\vec{\xi}$, $\vec{\eta}$). This frame is defined by its angle with respect to the one induced by the local flow motion that we denote as ($\vec{t}$, $\vec{n}$), with $\phi$ the angle of $\vec{t}$ with respect to the x-axis. This approach leads to simpler expressions for the mechanical equations. Some orientations of the vector $\vec{\xi}$ make the set of equations ill posed: $\vec{\xi} = \vec{t}$ for all flow regimes and angle($\vec{\xi}$, $\vec{t}$) = $\epsilon$ arcsin(1/M) $\epsilon = \pm 1$, if the flow is locally supersonic. The curves where at each point the tangent vector has one of these specific orientations with respect to the velocity are, respectively, the aforementioned streamtraces S, C$^+$ (for $\epsilon = 1$) and C$^-$ (for $\epsilon = -1$). Along these curves, equations with only $\vec{\xi}$ derivatives are satisfied. These ODEs are the classical DCEs.

For a 2D flow of an ideal gas (Section 2.1), without any assumptions of constant stagnation enthalpy, or constant entropy or null vorticity, (a) the stagnation enthalpy and
entropy are found to be constant along streamtraces; and (b) one DCE is found along the C$^+$C$^-$[3]. This DCE, expressed in terms of variations in $\phi$, $||\overline{U}||$, S, and H, reads

$$ \epsilon d\phi - \sqrt{M^2 - 1} \frac{d||\overline{U}||}{||\overline{U}||} - \sqrt{M^2 - 1} \frac{TdS - dH}{||\overline{U}||^2} = 0, \quad (14) $$

where the Crocco equation can be used to rewrite the last term. If H is constant over the whole fluid domain, the DCE also reads

$$ \epsilon d\phi - \frac{\sqrt{M^2 - 1}}{1 + \frac{\gamma - 1}{2} M^2} \frac{dM}{M} - \frac{\sqrt{M^2 - 1}}{\gamma r M^2} dS = 0. \quad (15) $$

If S is also constant, the flow is irrotational according to Crocco's theorem and the DCE reads

$$ \epsilon d\phi - \frac{\sqrt{M^2 - 1}}{1 + \frac{\gamma - 1}{2} M^2} \frac{dM}{M} = 0. \quad (16) $$

This equation can be integrated using the Prandtl-Meyer function,

$$ v(M) = \sqrt{\frac{\gamma + 1}{\gamma - 1}} \tan^{-1} \left( \sqrt{\frac{\gamma - 1}{\gamma + 1}} (M^2 - 1) \right) - \tan^{-1}(\sqrt{M^2 - 1}), $$

and the final algebraic equations resulting from the integration of (16) are

$$ k^{-} = \phi + v(M) \text{ is constant along a } C^{-} \qquad k^{+} = \phi - v(M) \text{ is constant along a } C^{+}. \quad (17) $$

Other sets of primitive variables may be used to derive (14) and the resulting simplified equations [4,31]. Also possible is the demonstration of these equations in a mapping of the plane rather than in the aforementioned frame of reference [31].

In case the flow is assumed to be irrotational from the beginning, a very fast demonstration is possible involving the velocity potential—see [32] or [2,34,35]. The governing nonlinear equation for a two-dimensional potential flow is:

$$ (1 - \frac{u^2}{c^2}) \frac{\partial u}{\partial x} - \frac{2uv}{c^2} \frac{\partial u}{\partial y} + (1 - \frac{v^2}{c^2}) \frac{\partial v}{\partial y} = 0 \quad (18) $$

and of course

$$ du = \frac{\partial u}{\partial x} dx + \frac{\partial u}{\partial y} dy \qquad dv = \frac{\partial v}{\partial x} dx + \frac{\partial v}{\partial y} dy. $$

This results in a simple 3 × 3 linear system for $(\partial u / \partial x)$, $(\partial v / \partial y)$, and $(\partial u / \partial y)$ (that is equal to $(\partial v / \partial x)$ for the potential flow) from which (16) and (17) are easily derived.

## 3.2. Derivation of the Characteristic Equations in Conservative Variables

The Cauchy problem aims at computing the derivatives $(\partial W / \partial x)$ and $(\partial W / \partial y)$ from the values of W at two neighboring points a and b. The relevant linear system has been posed in Section 2.2 in Equation (4). Its determinant is easily calculated as

$$ D = \begin{vmatrix} dxI & dyI & dxI & 0 \\ & & = & \\ A & B & A & B - (dy/dx)A \end{vmatrix} = dx^4 |B - (dy/dx)A| \quad (19) $$
Of course, $D = |dxB - dyA| = |-dxB + dyA|$ and the problem is ill posed along the S, C+, and C- curves, as recalled in Section 2.3. The fact that the Cauchy problem is ill posed along these specific curves, whereas $(\partial W/\partial x)$ and $(\partial W/\partial y)$ are actually defined and bounded, implies that not only the denominator in the Cramer formula applied to (4) is equal to zero, but also the eight numerators corresponding to the unknowns $(\partial\rho/\partial x)$, $(\partial\rho u/\partial x)$, $(\partial\rho v/\partial x)$, $(\partial\rho E/\partial x)$, $(\partial\rho/\partial y)$, $(\partial\rho u/\partial y)$, $(\partial\rho v/\partial y)$, and $(\partial\rho E/\partial y)$ must be equal to zero. This is precisely the principle of the method of characteristics that allows ODEs to be derived along the characteristic curves.

The principle of the calculation of the ODE satisfied along the S, C+, and C- curves is recalled for variable $(\partial\rho/\partial x)$, the boundedness of which requires that

$$ \begin{vmatrix} d\rho & 0 & 0 & 0 & dy & 0 & 0 & 0 \\ d\rho u & dx & 0 & 0 & dy & 0 & 0 & 0 \\ d\rho v & 0 & dx & 0 & 0 & dy & 0 & 0 \\ d\rho E & 0 & 0 & dx & 0 & 0 & dy & 0 \\ | & | & | & | & | & | & | & | \\ 0 & A_2 & A_3 & A_4 & B_1 & B_2 & B_3 & B_4 \end{vmatrix} = 0. \quad (20) $$

The determinant is expanded along the first column with the notation of the next equation:

$$ K_{1x}^1 d\rho - K_{1x}^2 d\rho u + K_{1x}^3 d\rho v - K_{1x}^4 d\rho E = 0. \quad (21) $$

In this equation, for example

$$ K_{1x}^4 = \begin{vmatrix} 0 & 0 & 0 & dy & 0 & 0 & 0 \\ dx & 0 & 0 & dy & 0 & 0 & 0 \\ 0 & dx & 0 & 0 & dy & 0 & 0 \\ | & | & | & | & | & | & | \\ A_2 & A_3 & A_4 & B_1 & B_2 & B_3 & B_4 \end{vmatrix} = \begin{vmatrix} 0 & 0 & 0 & dy & 0 & 0 & 0 \\ dx & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & dx & 0 & 0 & 0 & 0 & 0 \\ | & | & | & | & | & | & | \\ A_2 & A_3 & A_4 & B_1 & B_2 - tA_2 & B_3 - tA_3 & B_4 \end{vmatrix} \quad (22) $$

Finally,

$$ K_{1x}^4 = dx^2 dy \mid A_4 (B_2 - tA_2) (B_3 - tA_3) (B_4) \mid = dx^2 dy \kappa \gamma_1 (u + tv) \quad (23) $$

At this step, there is no assumption on the value of $t$ with respect to the velocity vector $(u, v)$. The other terms of the differential form are equal to

$$ K_{1x}^1 = -dx^3 \kappa (\kappa^2 v - \gamma_1 v (1+t^2) H + \gamma_1 (tu+v) E_c + 2 \gamma_1 t^2 v E_c) \quad (24) $$

$$ K_{1x}^2 = -dx^2 dy \kappa (\gamma u(u+tv) - 2E_c) \quad (25) $$

$$ K_{1x}^3 = dx^2 dy \kappa (\gamma v (u + t v) - 2 t E_c) \quad (26) $$

The characteristic Equation (21) can, thus, be rewritten as

$$ \begin{aligned} & - dx^3 \kappa (\kappa^2 v - \gamma_1 v (1 + t^2) H + \gamma_1 (tu + v) E_c + 2 \gamma_1 t^2 v E_c) d\rho \\ & + dx^2 dy \kappa (\gamma u(u + tv) - 2E_c) d\rho u + dx^2 dy \kappa (\gamma v (u + t v) - 2 t E_c) d\rho v \\ & - dx^2 dy \kappa \gamma_1 (u + tv) d\rho E = 0 \end{aligned} \quad (27) $$

If $dx \neq 0$ and $\kappa \neq 0$, this equation may be further simplified:

$$ \begin{aligned} & - (\kappa^2 v - \gamma_1 v (1 + t^2) H + \gamma_1 (tu + v) E_c + 2 \gamma_1 t^2 v E_c) d\rho \\ & + t (\gamma u(u + tv) - 2E_c) d\rho u + t (\gamma v (u + t v) - 2 t E_c) d\rho v \\ & - t \gamma_1 (u + tv) d\rho E = 0 \end{aligned} \quad (28) $$
Of course, along the S curves, $κ$ is null and Equation (28) seems to be only valid for the $C^+$ and $C^-$ curves. Actually, as $(-vdx + udy = κdx)$ has a multiplicity of two in the determinant D, (28) is also needed for the existence of $(∂ρ/∂x)$ along the S curves (see Section 3.3).

For the other seven partial derivatives, the expressions of the coefficients of the CEs are presented in Appendix A. In addition, it has been proved in Section 2 that, $W_l$ being one of the four conservative variables, the equations for the boundedness of $(∂W_l/∂x)$ and $(∂W_l/∂y)$ along the characteristics are proportional by a factor of $-t$. Therefore, only the counterparts of Equation (28) for the existence of $(∂ρu/∂x)$, $(∂ρv/∂x)$, and $(∂ρE/∂x)$ are presented hereafter:

$$ 
\begin{aligned}
& - \gamma_1 t E_c(-\gamma_1 H + (\gamma + 1) E_c) d\rho \\
& + (\gamma_1(\gamma t u + \kappa)E_c + v^2 \kappa + \gamma_1(v - \gamma_1 t u)H) d\rho u \\
& - t(u^2\kappa + \gamma_1(\kappa - \gamma v)E_c + \gamma_1(\gamma_1 v - u t)H) d\rho v \\
& + t \gamma_1(\gamma_1 H - (\gamma + 1)E_c) d\rho E = 0
\end{aligned} 
\quad (29) $$

$$ 
\begin{aligned}
& + \gamma_1 E_c t^2(\gamma_1 H - (\gamma + 1) E_c) d\rho \\
& - t(-v^2\kappa - \gamma_1\kappa E_c - \gamma\gamma_1 t u E_c + \gamma_1(\gamma_1 t u - v)H) d\rho u \\
& + (-v^3 + 3t u v^2 - 2t^2 u^2 v + \gamma_1((2 - \gamma)t^2 v - \kappa)H + \gamma_1(\gamma t^2 v + \kappa)E_c) d\rho v \\
& - t^2\gamma_1(-\gamma_1 H + (\gamma + 1)E_c) d\rho E = 0
\end{aligned} 
\quad (30) $$

$$ 
\begin{aligned}
& - \gamma_1 t(u + tv)E_c(\gamma_1 E_c + (2 - \gamma)H) d\rho \\
& + t(\gamma_1 u^2 + v\kappa + \gamma_1 t u v)(\gamma_1 E_c + (2 - \gamma)H) d\rho u \\
& - t(u\kappa - \gamma_1 t v^2 - \gamma_1 u v)(\gamma_1 E_c + (2 - \gamma)H) d\rho v \\
& + (-v\kappa^2 - \gamma_1((1 + \gamma t^2)v + \gamma_1 t u)E_c + \gamma_1^2 t(u + tv)H - \gamma_1\kappa H) d\rho E = 0
\end{aligned} 
\quad (31) $$

### 3.3. Ordinary Differential Equations along the Streamtraces S

Along the trajectories, $udy - vdx = 0$ is a zero of the denominator of Cramer's formulas with a multiplicity of two. The term $κdx = udy - vdx$ appears to be a factor of all $K_{mx}^l$ and the expressions obtained by removing this term from $K_{mx}^l$ are denoted by $\bar{K}_{mx}^l$. Obviously wherever the Cauchy problem is well defined

$$ \frac{\partial \rho}{\partial x} = \frac{\bar{K}_{1x}^1 d\rho - \bar{K}_{1x}^2 d\rho u + \bar{K}_{1x}^3 d\rho v - \bar{K}_{1x}^4 d\rho E}{(-v dx + u dy)(-v dx + u dy + c dl)(-v dx + u dy - c dl)}. \quad (32) $$

If the point a is fixed and the neighboring point b is moved closer and closer to $S_a$, $(-v dx + u dy) \to 0$ and the boundedness of $(∂ρ/∂x)$ requires the numerator of (32) to be equal to zero:

$$ \bar{K}_{1x}^1 d\rho - \bar{K}_{1x}^2 d\rho u + \bar{K}_{1x}^3 d\rho v - \bar{K}_{1x}^4 d\rho E = 0. \quad (33) $$

The expressions for the $\bar{K}$ coefficients are much simpler than those of the $K$ coefficients thanks to the nullity of $κ$. It can be verified by hand or by formal calculation that the CEs expressing the boundedness of $(∂W_l/∂x)$ and their counterparts for $(∂W_l/∂y)$ are proportional by a factor of $-t$. It is also possible to check that $\bar{K}_{1x}^l = -t\bar{K}_{ly}^l$ and then use the continuity of the $\bar{K}$ coefficients as functions of $(W, dx, dy, \gamma)$ and the validity of $\bar{K}_{mx}^l = -t\bar{K}_{my}^l$ $l \gg m$ outside the variety $udy - vdx = 0$. We denote $(\gamma_1 E_c + (2 - \gamma)H)$ as G; the fully simplified differential forms then read

$$ (2E_c - H) d\rho - u d\rho u - v d\rho v + d\rho E = 0 \quad (34) $$

$$ ((\gamma + 1)E_c - \gamma_1 H) (E_c d\rho + d\rho E) - (\gamma E_c + (2 - \gamma)H) (ud\rho u + v d\rho v) = 0 \quad (35) $$

$$ -E_c G d\rho + uG d\rho u + vG d\rho v + (\gamma_1 H - \gamma E_c) d\rho E = 0, \quad (36) $$
where Equation (35) is obtained twice, from the existence of the derivatives of both momentum components. The rank of the matrix was first calculated with the Maple software 2020.0 (Supplementary Materials) and, as expected, is two. It can also be verified by hand calculation that $Ec \times (34)$ minus (35) is exactly Equation (36) and then, finally, that the $d\rho$ and $d\rho E$ columns in {(34),(35)} are not linked.

Along the streamtraces, the derivative of the total enthalpy, $H$, and the entropy, $S$, are expected to be null [27,32]. Before verifying that these properties are ensured by the CE found just before, the differential of $p/\rho^\gamma$ times $\rho^\gamma$ and the differential of $H$ times $\rho$ are calculated:

$$ \rho^{\gamma} d(p/\rho^{\gamma}) = (2E_c - H) d\rho - u d\rho u - v d\rho v + d\rho E \quad (37) $$

$$ \rho dH = (-\gamma E + 2(\gamma - 1)E_c) d\rho - \gamma_1 u d\rho u - \gamma_1 v d\rho v + \gamma d\rho E. \quad (38) $$

Equation (34) is exactly the same as $\rho^\gamma d(p/\rho^\gamma) = 0$. This yields $d(p/\rho^\gamma) = 0$ and obviously the derived CEs contain the property of constant entropy along the streamtraces. In addition, $G \times (34)$ plus (36) yields, after simplification by a common factor, $-Gd\rho + d\rho E = 0$. Finally, adding this equation to $\gamma_1$ times Equation (34), precisely yields $\rho dH = 0$, so that the second expected property for Euler flows along the trajectories is contained in the derived CEs. Of course, the conservation of the stagnation enthalpy and the entropy along the trajectories of a steady-state inviscid flow is also verified in 3D [27,32].

### 3.4. Ordinary Differential Equations along the $C^+C^-$

The proportionality of all eight differential forms along the $C^+$ and $C^-$ has been proved in Section 2.4 and Equation (28) may be retained as the relevant DCE. The factor of $d\rho$ is quite complex in (28) but it admits a simpler form derived from the degree two equations satisfied by $t^+$ and $t^-$: these curve slopes are defined as

$$ t^{\pm} = \tan(\phi + \beta) \quad \text{with} \quad \tan \phi = \frac{v}{u} \quad \sin \beta = \pm \frac{1}{M}, \quad (39) $$

or alternatively by their explicit expressions [30]:

$$ t^{\pm} = \frac{uv \pm c\sqrt{u^2 + v^2 - c^2}}{u^2 - c^2}. \quad (40) $$

They are the two roots of the following degree two equation [30]

$$ \gamma_1 (1 + t^2) H = \gamma_1 (1 + t^2) E_c + (tu - v)^2, \quad (41) $$

which allows the expression of (28) to be simplified along a $C^+$ or a $C^-$. The final DCE reads

$$ \gamma_1(u + tv)E_c d\rho - (\gamma u(u + tv) - 2E_c) d\rho u - (\gamma v(u + tv) - 2tE_c) d\rho v + \gamma_1(u + tv) d\rho E = 0 \quad (42) $$

with $t = t^+$ for a $C^+$ and by $t = t^-$ for a $C^-$.

It has been first checked that this $C^+$ (resp. $C^-$) differential form is linked with $\{dk^+, dH, dS\}$ (resp. $\{dk^-, dH, dS\}$) consistently with Equation (14). It is proved in Appendix B that these CEs, expressed with different variables, are the same.

## 4. Assessment of the Direct and Adjoint Characteristic Equations for a Nozzle Flow

### 4.1. Supersonic Nozzle Configuration

The test case for the assessment of the DCEs and ACEs is a supersonic nozzle designed for an aircraft flying at Mach 1.6. The geometry was defined in the framework of the SENECA EU project [36]. The nozzle is a classical convergent-divergent duct. At the subsonic injection plane, the stagnation conditions are fixed as well as the direction of the velocity (aligned with the x-axis). The ratios of the far-field flow total pressure (resp.
temperature) to the total pressure (resp. total temperature) at the inlet are $p_{i_{inl}}/p_{i_{\infty}} = 1.605$ and $T_{i_{inl}}/T_{i_{\infty}} = 1.606$. As an inviscid flow is calculated, these two conditions and the far-field Mach number fully define the flow. Figure 1 shows a structural diagram of the case object and the fluid domain.

<figure><img src="image_4.png" /><figcaption>Figure 1. SENECA nozzle: fluid domain and boundary conditions diagram.</figcaption></figure>

The original geometry was designed for a circular axisymmetric engine and the nozzle contour in the symmetry plane was simply extracted to define the 2D geometry. Consequently, the mass flow rates are different. Due the different laws of areas along the x-coordinate (the inlet and main flow direction), the flows are also different along the midline of the nozzle. Nevertheless, the Mach number distribution along this line is similar, with $M = 1$ at the throat and $M \simeq 1.9$ at the exit.

The functional output of interest is the thrust $Th$,

$$Th = \int_{\Gamma_w} (p - p_{\infty})ds_x + \int_{\Gamma_{inl}} (\rho u^2 + p)ds_x, \quad (43)$$

where $\Gamma_w$ is the internal wall of the nozzle, $\Gamma_{inl}$ the injection section—see Figure 2—, and $s_x$ the x-component of the surface vector (oriented towards the outside of the fluid domain).

<figure><img src="image_5.png" /><figcaption>Figure 2. SENECA nozzle. Green curve: $C^+$ passing by the trailing edge of the nozzle. Dark gray zone: zone of the flow influencing the thrust $Th$.</figcaption></figure>

## 4.2. Flow Calculation, Flow Analysis

A structured eight-block mesh with 370,538 cells discretizes the fluid domain. The flow simulations are run with the elsA code [37] using the Jameson–Schmidt–Turkel (JST) scheme [38].

elsA is a finite-volume cell-centered code that generally requires specific state variable calculations at the boundary faces. The exception is the symmetry line, where the centered and artificial dissipation fluxes are calculated with current-face formulas involving ghost cells across the symmetry line. At the subsonic inlet and at the non-reflecting boundaries, no artificial dissipation flux is added. The physical inviscid flux is evaluated there for state variables that are most often derived from the classical theory of characteristics—see chapter 3 and chapter 11 in [27]. At the non-reflecting boundaries, its equations are applied using switches based on the sign of the eigenvalues of the normal Jacobian. At the subsonic inlet boundary, one characteristic equation corresponds to convection from
the fluid domain towards the boundary and two numerical approaches are possible with elsA: (a) the first one, based on characteristics, calculates the static pressure $p_b$ and the velocity $u_b$ at the inlet boundary by equalizing ($p - \bar{p}\bar{c}u$) between the boundary and the fluid domain. In addition, $p_b$ and $u_b$ are linked by a second equation derived from the fixed stagnation state—$p_b = p_i(1 - u_b^2/(2c_pT_i))^{\frac{\gamma}{\gamma-1}}$— which results in a nonlinear equation in $u_b$ that is solved using a Newton method; (b) the second one extrapolates the normal velocity from the fluid domain to the inlet boundary and forces a value in the interval of subsonic velocities compatible with the fixed total enthalpy. All other variables at the inlet boundary are then derived from the stagnation state. The latter approach is more robust for nozzle flows and is the one selected for the current calculation.

The accuracy of the discrete flow is estimated through an $L_1$ relative error calculation for the total enthalpy, for which exact theoretical constant values are known for both the jet plume and the external flow. This error is found to be $2.3 \times 10^{-3}$ for the jet and $2.5 \times 10^{-3}$ for the external flow.

For detailed flow analysis, Figure 3 shows the Mach number versus static pressure isolines in the computational and mirrored domain, with $C^+$ and $C^-$ emanating from the nozzle trailing edges and also the highest Mach number isolines in a specific color (the full length of the domain in $x$ appears in Figures 4 and 5; the full length of the domain in $y$ is mentioned in Figure 1).

A slip line emanating from the trailing edge of the nozzle (which corresponds to an S) allows the jet plume shape to be clearly defined—see Figure 4. For this inviscid flow, there is no significant shear layer growth as the fluid goes downstream of the nozzle and the numerical shear layer becomes thinner and thinner as the mesh is refined. It is then possible to post-process any flow data inside or outside the jet plume, such as the mass flow rate, the entropy losses, the thrust accounting, or any data relying on the integration of values inside the jet.

Outside the nozzle, a shock attached to the nozzle trailing edge follows a $C^+$ on the top side (and a $C^-$ on the opposite side). That shock compresses the flow and increases the static pressure to approximately 14,000 Pa in the vicinity of the nozzle’s trailing edge. Then, the flow extends again, and a second outside shock, much weaker, is located at $x \approx 2.7$ m. This is due to the interaction between the two expansion fans inside the nozzle that intersect at $x \approx 2.0$ m, generating a shock in the inner flow (clearly visible in the static pressure figure). That shock crosses the jet plume. Finally, the flow accelerates again in the expansion fan up to the far-field static pressure, around 9700 Pa.

Inside the nozzle, when exiting at $x \approx 1.2$ m, the jet is supersonic and over-pressurized, so it extends into a lower pressure with a diverging shape. Being divergent and supersonic, the jet Mach number increases, up to the point at $x \approx 2.7$ m where the aforementioned shock compresses the jet. Then, the expansion fan from the nozzle trailing edge expands the flow once again, up to a pressure close to the far-field static pressure. The Mach number of the jet is then nearly constant, close to $M = 1.9$. Some isolines with a dense scale are plotted in white in Figure 3 in order to further highlight the weak interactions for $x > 2.7$ m. Four diamonds are limited by $C^+$ and $C^-$. In these zones, the flow is slightly perturbed (also visible in the static pressure figure). Therefore, the theory of the characteristic equations is suitable to analyze this Euler supersonic flow.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** SENECA nozzle. Isolines of Mach number (with Mach numbers in [1.9, 2.0] in white) vs. isolines of static pressure. Including selected $C^+$ in green and $C^-$ in purple for detailed flow analysis.</figcaption></figure>
### 4.3. Thrust-Adjoint Calculation: Extraction of C⁻ and C⁺ Curves

The adjoint calculation for the thrust is run with the discrete adjoint module of the elsA code [39], which is used for both research [12,40] and application purposes. In [40], it was proved that the discrete adjoint of the JST scheme on a structured mesh is dual consistent for Euler flows, but in the cells next to those adjacent to a boundary. This adjoint consistency property inside the fluid domain is very valuable to numerically assess continuous adjoint properties. The change in the Jacobian to obtain the dual consistency in the vicinity of the physical boundaries is complex to implement in an industrial code. In this study, it is not used and slight oscillations in the adjoint field are observed near the boundaries. For the sake of simplicity, considering the mesh density, the points inside the first two cells in the vicinity of a boundary are removed from the extracted characteristic curves. These curves were extracted from a shared point and the field of tangent vectors defined by (1, t)—more precisely, (1, v/u) for the S, (1, t⁺) for the C⁺, and (1, t⁻) for the C⁻. The common point is located in (615,130), in the middle of the x-range of the divergent duct and half of the height of the upper half of the channel at the abscissa—see Figure 4. Sufficiently large numbers of points have been extracted along the characteristic curves to obtain a converged error analysis, presented in Sections 4.4 and 4.5.

The first component of the adjoint field is presented in Figure 5. The exact adjoint of the thrust is expected to be zero downstream of the backward C⁺ characteristic (and corresponding to C⁻ in the symmetric non-plotted domain) emanating from the rear of the nozzle, as no perturbation downstream of these lines can affect the flow on the support of T (we recall that the velocity is supersonic in the divergent part of the nozzle). This property is well satisfied by the discrete adjoint.

<figure><img src="image_7.png" /><figcaption>Figure 4. SENECA nozzle. Isolines of Mach number. Selected S (in black), C⁺ (in green), and C⁻ (in purple) for the validation of the DCEs and ACEs. The thicker S (in dark gray) delimits the internal and the external flow.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 5. SENECA nozzle. Isolines of the discrete adjoint of the thrust, component $\psi_1$ (dual of the mass-flow residual).</figcaption></figure>

### 4.4. Relation between the Adjoint Components of Density Equation and Total Energy Equation

One of the few known properties of the Euler exact adjoint fields—denoted here as $\lambda$—is that along a trajectory

$$ \frac{d\lambda_1}{ds} - H\frac{d\lambda_4}{ds} = 0. \qquad (44) $$

This equation can very simply be derived from the two streamtraces of the ACEs,

$$ \frac{d\lambda_1}{ds} + u \frac{d\lambda_2}{ds} + v \frac{d\lambda_3}{ds} + E_c \frac{d\lambda_4}{ds} = 0 \qquad E_c \frac{d\lambda_1}{ds} + H(u \frac{d\lambda_2}{ds} + v \frac{d\lambda_3}{ds}) + H^2 \frac{d\lambda_4}{ds} = 0 \quad (45) $$

or can be directly proved by more complex calculations from the Euler adjoint continuous equations [30]. As the stagnation enthalpy H is constant along a trajectory, Equation (44) yields $\lambda_1 = H\lambda_4 + \lambda^*$ with, a priori, a constant $\lambda^*$ depending on the streamtrace. However, Giles and Pierce proved, using a physical term source approach, that the simpler formula

$$ \lambda_1 = H\lambda_4 \qquad (46) $$
holds in the common case where the function of interest depends only on the static pressure [41]. For our nozzle application, all trajectories, when traveled in the adjoint information propagation direction (that is, backwards) start in a zone where the adjoint vector is zero, so that (46) is true over the whole fluid domain, although the nozzle thrust does not exclusively depend on the static pressure. We expect the corresponding equation at the discrete level, $\psi_1 = H\psi_4$, to be well satisfied; this is indeed the case—see Figure 6.

A series of recent publications have studied the properties the Euler lift- and drag-adjoint fields, in particular in the vicinity of walls [40,42–45]. How accurately Equation (46) is satisfied at a discrete level is discussed in [40] and, more generally, the empirical point of view on code verification as put forward in [33] §3.2.1 seems very relevant to us, and we see Equations (45) and (46) as both theoretical results and useful tools for code verification.

<figure><img src="image_9.png" /><figcaption>Figure 6. SENECA nozzle. Isolines of $\psi_1$ (left) and $H \times \psi_4$ (right) with same levels.</figcaption></figure>

### 4.5. Assessment of the DCEs

The numerical assessment method consists in integrating the DCEs along the characteristic curves for the aforementioned fine-mesh flow. More precisely, for the $C^+$ and $C^-$, the following integrals are calculated using (42)

$$ \Xi C^+ = \int_{C^+} \left( \gamma_1(u + t^+ v) E_c \frac{d\rho}{ds} - (\gamma u(u + t^+ v) - 2E_c) \frac{d\rho u}{ds} - (\gamma v(u + t^+ v) - 2t^+ E_c) \frac{d\rho v}{ds} + \gamma_1(u + t^+ v) \frac{d\rho E}{ds} \right) ds $$

$$ \Xi C^- = \int_{C^-} \left( \gamma_1(u + t^- v) E_c \frac{d\rho}{ds} - (\gamma u(u + t^- v) - 2E_c) \frac{d\rho u}{ds} - (\gamma v(u + t^- v) - 2t^- E_c) \frac{d\rho v}{ds} + \gamma_1(u + t^- v) \frac{d\rho E}{ds} \right) ds, $$

with s the curvilinear abscissa along the characteristic curve. The subparts of $\Xi C^+$ and $\Xi C^-$ are also calculated to avoid any error in scale when discussing close to zero numerical values. For $\Xi C^+$, for example,

$$ \Xi C_1^+ = \int_{C^+} \left( \gamma_1 E_c (u + vt^+) \frac{d\rho}{ds} \right) ds \qquad \Xi C_2^+ = - \int_{C^+} \left( (\gamma u(u + vt^+) - 2E_c) \frac{d\rho u}{ds} \right) ds $$

$$ \Xi C_3^+ = - \int_{C^+} \left( (\gamma v(u + vt^+) - 2E_c t^+) \frac{d\rho v}{ds} \right) ds \qquad \Xi C_4^+ = \int_{C^+} \left( \gamma_1(u + vt^+) \frac{d\rho E}{ds} \right) ds. $$

Of course, the $\Xi C^+ = \Xi C_1^+ + \Xi C_2^+ + \Xi C_3^+ + \Xi C_4^+$ (resp. $\Xi C^- = \Xi C_1^- + \Xi C_2^- + \Xi C_3^- + \Xi C_4^-$) sum is expected to be much smaller than its subparts. The integration is performed in the information propagation direction (increasing x) along the selected $C^+$ and the selected $C^-$ presented in Figure 4. Almost null values of $\Xi C^+$ and $\Xi C^-$ are indeed observed—see Figure 7.

<figure><img src="image_10.png" /><figcaption>Figure 7. Numerical assessment of DCE (42). $\Xi C^+$ (left), $\Xi C^-$ (right) and their subparts integrated along the selected $C^+ C^-$. Method of verification: the black curve should ideally coincide with the x-axis.</figcaption></figure>
An error per series of curves has been calculated by computing the max of $|\Xi C^+|$ (resp. $|\Xi C^-|$) divided by the max of the absolute value of all four corresponding subparts $|\Xi C_l^+|_{l \in \{1,2,3,4\}}$ (resp. $|\Xi C_l^-|_{l \in \{1,2,3,4\}}$). The errors are $9.0 \times 10^{-8}$ for $\Xi C^+$ and $4.2 \times 10^{-8}$ for $\Xi C^-$. Equivalent successful verifications have been performed to assess the validity of (34) and (35). They are not presented here for the sake of brevity. These two streamtrace DCEs are equivalent to $dS = 0$ and $dH = 0$, as discussed in Section 3.3. As a complementary verification, the relative variation in $H$ (over the whole streamtrace) and $S$ (up to the first shockwave) are calculated. Relative errors of $2.3 \times 10^{-6}$ and $1.4 \times 10^{-6}$ are found.

## 4.6. Assessment of the ACEs for the Thrust Adjoint

For the thrust adjoint, the assessment method is the same except that the ACEs are integrated in the forward sense for the adjoint, that is, from the support of the functional output backwards with respect to the direction of the flow information propagation. The integrals to evaluate (as well as their subparts) from the discrete flow and thrust-adjoint fields are [30]:

$$ \Xi_{adj}S^1 = \int_S \left( E_c \frac{d\psi_1}{ds} + H(u \frac{d\psi_2}{ds} + v \frac{d\psi_3}{ds}) + H^2 \frac{d\psi_4}{ds} \right) ds \quad (47) $$

$$ \Xi_{adj}S^2 = \int_S \left( \frac{d\psi_1}{ds} + u \frac{d\psi_2}{ds} + v \frac{d\psi_3}{ds} + E_c \frac{d\psi_4}{ds} \right) ds \quad (48) $$

$$ \Xi_{adj}C^+ = \int_{C^+} \left( (u + vt^+) \frac{d\psi_1}{ds} + (u^2 + v^2) \left( \frac{d\psi_2}{ds} + t^+ \frac{d\psi_3}{ds} \right) + H(u + vt^+) \frac{d\psi_4}{ds} \right) ds \quad (49) $$

$$ \Xi_{adj}C^- = \int_{C^-} \left( (u + vt^-) \frac{d\psi_1}{ds} + (u^2 + v^2) \left( \frac{d\psi_2}{ds} + t^- \frac{d\psi_3}{ds} \right) + H(u + vt^-) \frac{d\psi_4}{ds} \right) ds. \quad (50) $$

As in Section 4.3, the CEs are numerically assessed by checking that $\Xi_{adj}S^1$, $\Xi_{adj}S^2$, $\Xi_{adj}C^+$ and $\Xi_{adj}C^-$ are much smaller than their subterms all along the characteristic curves. This property is actually well satisfied, as can be seen in Figure 8. The error per series of curves, as defined in the previous subsection, has been calculated: $1.8 \times 10^{-6}$ for $\Xi_{adj}S^1$, $2.6 \times 10^{-6}$ for $\Xi_{adj}S^2$, $2.0 \times 10^{-6}$ for $\Xi_{adj}C^+$, and $9.5 \times 10^{-9}$ for $\Xi_{adj}C^-$. Method of verification: the black curve should ideally coincide with the x-axis.

<figure><img src="image_11.png" /><figcaption>Figure 8. Numerical assessment of the ACEs [30]. $\Xi_{adj}S^1$ (upper left), $\Xi_{adj}S^2$ (upper right), $\Xi_{adj}C^+$ (bottom left), $\Xi_{adj}C^-$ (bottom right) and their subparts integrated along the selected $S$, $C^+$, and $C^-$. Method of verification: the black curve should ideally coincide with the x-axis.</figcaption></figure>
As mentioned in [30], a lower accuracy would be observed for more upwind C⁺ and C⁻. Our understanding of this stiffness is the sensitivity of the $c\sqrt{u^2 + v^2 - c^2} = c^2\sqrt{M^2 - 1}$ term in the slope of Equation (40) close to the boundary of the supersonic zone, $M = 1^+$.

## 5. Conclusions

The well-known equations satisfied in 2D inviscid flows along trajectories and in supersonic zones along C⁺ and C⁻, are very helpful for the understanding of Mach waves, shockwaves, rarefaction waves, and nozzle flows. In this paper it has been shown that these equations can be derived using the conservative variables in a Cartesian frame of reference. This is in contrast to the canonical demonstrations in which the primitive variables are used in a frame of reference relative to the axis of motion.

In this work the study of the characteristic equations of the adjoint Euler equations has been continued, in particular using linear algebra arguments that simplify the discussion of the number of independent equations satisfied along the C⁺ and C⁻, also deriving consequences of the characteristic equations relative to the adjoint components of the density and total energy equations. Of course, most of the efforts in CFD are devoted today to RANS, LES, and DNS flows. However, exact equations for Euler adjoint fields satisfied along trajectories along C⁺ and C⁻ in supersonic zones or throughout the fluid domain are very valuable for code verification. In the case of a continuous adjoint scheme or dual consistent discrete adjoint, they may directly serve for code verification and even validation of an adjoint code module. In the case of a general discrete adjoint code, the assessment of numerical solutions with respect to exact continuous equations is also linked to the difficult question of adjoint consistency, about which little has been proved today.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10.3390/aerospace10090797/s1, three Python files checking the formulas of $K_{mx}^l, K_{my}^l, \bar{K}_{mx}^l, \bar{K}_{my}^l$. Two Maple files checking the rank of the differential forms along the S and C⁺ curves.

**Author Contributions:** Investigation, K.A. and J.P.; software, J.P.; validation, K.A. and O.A.; writing, J.P., K.A. and O.A.; formal analysis, K.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** The nozzle considered in Section 4 was designed in the framework of the SENECA project that is funded by the European Union's Horizon 2020 research and innovation program under grant agreement No. 101006742, project SENECA (LTO) Noise and Emissions of Supersonic Aircraft). This work was partially supported by the SONICE project, granted by the French Directorate General for Civil Aviation (DGAC).

**Data Availability Statement:** The data that support the findings of this study may be shared. Each request will be examined in conjunction with the managers of SENECA projects.

**Acknowledgments:** elsA V5.0 used in Section 4 and Appendix C is ONERA-Safran property. The authors warmly thank Sébastien Heib, Sana Amri, Sébastien Bourasseau, Thomas Hennion, and Antoine Dumont for many fruitful discussions.

**Conflicts of Interest:** The authors declare no conflict of interest.

## Appendix A

* The $K_{2x}^l$ coefficients are expressed below

$$K_{2x}^1 = -dx^3 \kappa \gamma_1 t E_c (\gamma_1 H - (\gamma + 1)E_c)$$

$$K_{2x}^2 = dx^3 \kappa (\gamma_1^2 tuE_c + \gamma_1 E_c \kappa + v^2 \kappa + \gamma_1 H v - \gamma_1^2 Htu + \gamma_1 tuE_c)$$

$$K_{2x}^3 = -dx^3 \kappa t (-u^2 \kappa + \gamma_1 ((\gamma + 1) v - t u) E_c + \gamma_1 H (t u - \gamma_1 v))$$

$$K_{2x}^4 = -dx^3 \kappa \gamma_1 t (-\gamma_1 H + (\gamma + 1)E_c)$$
• The $K_{3x}^l$ coefficients are expressed below

$$K_{3x}^1 = -dx^3 \kappa \gamma_1 t^2 E_c (-\gamma_1 H + (\gamma+1)E_c)$$

$$K_{3x}^2 = -dx^3 \kappa t (\gamma_1 ((\gamma+1) t u - v) E_c + v^2 \kappa + \gamma_1 H (v - \gamma_1 t u))$$

$$K_{3x}^3 = dx^3 \kappa (-v^3 + 3 t u v^2 - \gamma_1 H \kappa - 2 t^2 u^2 v + \gamma_1^2 t^2 v E_c + \gamma_1 H t^2 v \\ - \gamma_1^2 H t^2 v + \gamma_1 t^2 v E_c + \gamma_1 E_c \kappa)$$

$$K_{3x}^4 = -dx^3 \kappa \gamma_1 t^2 (\gamma_1 H - (\gamma+1) E_c)$$

• The $K_{4x}^l$ coefficients are expressed below

$$K_{4x}^1 = dx^3 \kappa \gamma_1 t E_c (u + t v) (\gamma_1 E_c + (2 - \gamma) H)$$

$$K_{4x}^2 = dx^3 \kappa t (2 \gamma_1 E_c + \gamma \kappa v) (\gamma_1 E_c + (2 - \gamma) H)$$

$$K_{4x}^3 = -dx^3 \kappa t (-u \kappa + \gamma_1 t v^2 + \gamma_1 u v) (\gamma_1 E_c + (2 - \gamma) H)$$

$$K_{4x}^4 = dx^3 \kappa (-\gamma_1 H \kappa - \gamma_1 v E_c - v \kappa^2 - \gamma_1^2 t u E_c - \gamma_1 t^2 v E_c \\ - \gamma_1^2 t^2 v E_c + \gamma_1^2 H t u + \gamma_1^2 H t^2 v)$$

• The $K_{1y}^l$ coefficients are expressed below

$$K_{1y}^1 = dx^3 \kappa (u \kappa^2 + 2 \gamma_1 u E_c + \gamma_1 t v E_c + \gamma_1 t^2 u E_c - \gamma_1 H (1 + t^2) u)$$

$$K_{1y}^2 = dx^3 \kappa (\gamma u(u + tv) - 2 E_c)$$

$$K_{1y}^3 = -dx^3 \kappa (\gamma v(u + tv) - 2tE_c)$$

$$K_{1y}^4 = -dx^3 \kappa \gamma_1 (t v + u)$$

• The $K_{2y}^l$ coefficients are expressed below

$$K_{2y}^1 = dx^3 \kappa \gamma_1 E_c (\gamma_1 H - (\gamma+1) E_c)$$

$$K_{2y}^2 = dx^3 \kappa (\gamma_1 E_c (t \kappa - \gamma u) + t^2 u^3 + 2uv^2 - 3 t u^2 v + \gamma_1 H ((\gamma - 2)u - t \kappa))$$

$$K_{2y}^3 = dx^3 \kappa (-u^2 \kappa + \gamma_1 ((\gamma+1) v - t u) E_c + \gamma_1 H (t u - \gamma_1 v))$$

$$K_{2y}^4 = dx^3 \kappa \gamma_1 (-\gamma_1 H + (\gamma+1) E_c)$$

• The $K_{3y}^l$ coefficients are expressed below

$$K_{3y}^1 = dx^3 \kappa \gamma_1 E_c t (-\gamma_1 H + (\gamma+1) E_c)$$

$$K_{3y}^2 = dx^3 \kappa (\gamma_1 ((\gamma+1) t u - v) E_c + v^2 \kappa + \gamma_1 H (v - \gamma_1 t u))$$

$$K_{3y}^3 = dx^3 \kappa t (\gamma_1 H (\gamma_1 v - t u) - 2 \gamma_1 v E_c - \gamma_1^2 v E_c + u^2 \kappa + \gamma_1 t u E_c)$$

$$K_{3y}^4 = dx^3 \kappa \gamma_1 t (\gamma_1 H - (\gamma+1) E_c)$$

• The $K_{4y}^l$ coefficients are expressed below

$$K_{4y}^1 = -dx^3 \kappa \gamma_1 E_c (u + t v) (\gamma_1 E_c + (2 - \gamma) H)$$

$$K_{4y}^2 = -dx^3 \kappa (2 \gamma_1 E_c + \gamma \kappa v) (\gamma_1 E_c + (2 - \gamma) H)$$

$$K_{4y}^3 = dx^3 \kappa (-u \kappa + \gamma_1 t v^2 + \gamma_1 u v) (\gamma_1 E_c + (2 - \gamma) H)$$

$$K_{4y}^4 = dx^3 \kappa (u \kappa^2 + \gamma \gamma_1 u E_c + \gamma_1 t^2 u E_c + \gamma_1^2 t v E_c - \gamma_1^2 H u - \gamma_1^2 H t v \\ - \gamma_1 H t^2 u + \gamma_1 H t v)$$
## Appendix B

Let us prove the equivalence between (14) and (42) in the case of a C⁺. Equation (14) reading

$$d\phi - \sqrt{M^2 - 1} \frac{d||\overline{U}||}{||\overline{U}||} - \sqrt{M^2 - 1} \frac{TdS - dH}{||\overline{U}||^2} = 0,$$

with

$$\phi = \arctan\left(\frac{\rho v}{\rho u}\right) \quad \text{and} \quad T = \frac{p}{\rho r} \quad \text{so that} \quad TdS = \frac{p}{\rho\gamma_1}d(\ln(\frac{p}{\rho\gamma})) = \frac{p}{\rho\gamma_1}(\frac{dp}{p} - \gamma\frac{d\rho}{\rho}).$$

A series of simplifications leads to the shorter expression

$$u \, d\rho v - v \, d\rho u + \sqrt{M^2 - 1} \, dp = 0.$$

Expressing $p$ in primitive variables yields

$$\gamma_1 E_c \sqrt{M^2 - 1} \, d\rho - (v + \gamma_1 u \sqrt{M^2 - 1}) \, d\rho u + (u - \gamma_1 v \sqrt{M^2 - 1}) \, d\rho v + \gamma_1 \sqrt{M^2 - 1} \, d\rho E = 0. \quad (A1)$$

In addition, Equation (42) reads

$$\gamma_1 E_c (u + t^+ v) \, d\rho - (\gamma u (u + t^+ v) - 2E_c) \, d\rho u - (\gamma v (u + t^+ v) - 2E_c t^+) \, d\rho v + \gamma_1 (u + t^+ v) \, d\rho E = 0 \quad (A2)$$

It appears that the coefficients of Equations (A1) and (A2) are proportional:

$$\frac{\gamma_1 E_c (u + t^+ v)}{\gamma_1 E_c \sqrt{M^2 - 1}} = \frac{-(\gamma u (u + t^+ v) - 2E_c)}{-(v + \gamma_1 u \sqrt{M^2 - 1})} = \frac{-(\gamma v (u + t^+ v) - 2E_c t^+)}{(u - \gamma_1 v \sqrt{M^2 - 1})} = \frac{\gamma_1 (u + t^+ v)}{\gamma_1 \sqrt{M^2 - 1}} = (u + t^+ v) \tan \beta. \quad (A3)$$

## Appendix C. Assessment of the Thrust Adjoint via a Shape Optimization

The discrete adjoint of the thrust is assessed via a shape optimization. In order to improve the thrust of the nozzle, six shape parameters are defined as factors of the following functions that drive changes in the y-coordinates:

$$D_y^l(x, y) = \frac{4(x - x_{l-1})^2(x_{l+1} - x)^2(y_H - (y_{wall}(x) - y))^2}{(x_{l+1} - x_{l-1})^4 y_H^2} \quad \text{for } x \in [x_{l-1}, x_{l+1}]$$

$$D_y^l(x, y) = 0 \quad \text{if } x \notin [x_{l-1}, x_{l+1}],$$

$l \in \{1, 6\}$. The seven points defining the x-range of the six bumps are linearly distributed between the abscissa of the inlet ($x_0$) and the abscissa of the rear of the nozzle ($x_7$). $y_H$ is equal to half of the height between the inside wall at the throat and the symmetry line. These $D_y^l$ functions allow the y-coordinates of the deformed meshes to be calculated from the design parameters $\alpha$ and the coordinates $x^0, y^0$ in the original mesh. The formula

$$y(\alpha_1 \dots \alpha_6) = y^0 + \alpha_1 D_y^1(x^0, y^0) + \alpha_2 D_y^2(x^0, y^0) + \alpha_3 D_y^3(x^0, y^0) + \alpha_4 D_y^4(x^0, y^0) \quad (A4)$$

$$+ \alpha_5 D_y^5(x^0, y^0) + \alpha_6 D_y^6(x^0, y^0), \quad (A5)$$

is applied point-wise and the x-coordinates in the deformed mesh are the same as in the original mesh.

In order to optimize the thrust, a steepest descent has been performed with a maximum allowed wall displacement of 5% of the throat height at each step. The norm of the first adjoint gradient is $||\nabla Th|| = 12.5 e^3$ with strong negative $\alpha_3$ and $\alpha_4$ components (this negative sign is because the thrust is negative when computed along the x-axis). The 1D-descent step is bounded by the geometric constraint. The resulting deformation is a significant enlargement of the throat—see Figures A1 and A3. This first optimization step increases the thrust by 3.1%.
<figure><img src="image_12.png" /><figcaption>Figure A1. SENECA nozzle. Isolines of the Mach number for the deformed geometry at step 1.</figcaption></figure>

A second optimization step is performed. The norm of the thrust gradient is $||\nabla Th|| = 12.3 e^3$. The displacement at this step is fixed by the 1D search of an optimum in the direction of the adjoint gradient. The nozzle is slightly enlarged in the inlet part—see Figures A2 and A3. The thrust is increased by 0.7%.

<figure><img src="image_13.png" /><figcaption>Figure A2. SENECA nozzle. Isolines of the Mach number for the deformed geometry at step 2.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure A3. SENECA nozzle. Inside wall shapes at the different steps: red, initial geometry; blue, after the first optimization step; black, after the second optimization step.</figcaption></figure>

Finally, we check the consistency of our results with the classical formulas for a quasi-1D convergent-divergent “de Laval” nozzle. The thrust for such a configuration is expressed as

$$Th_{\overline{1d}} = \dot{m}_{\overline{1d}} U_{exit} = \dot{m}_{\overline{1d}} M_{exit} \sqrt{\gamma r T_{exit}} \\ Th_{\overline{1d}} = \dot{m}_{\overline{1d}} M_{exit} \sqrt{\frac{\gamma r T_{inl}}{1 + 0.5 \gamma_1 M_{exit}^2}} = \dot{m}_{\overline{1d}} \sqrt{\frac{\gamma r T_{inl}}{1/M_{exit}^2 + 0.5 \gamma_1}},$$

where $\dot{m}_{\overline{1d}}$ is the mass flow. With a fixed stagnation state at the inlet and a value of $M_{exit}$ significantly higher than 1, we expect the thrust to be optimized only via an increase of the mass flow.

In addition, the quasi-1D nozzle flows are calculated from the well-known equation expressing the mass flow as a function of the local section Y and the corresponding Mach number:

$$\dot{m}_{\overline{1d}} = \frac{p_{inl}}{\sqrt{T_{inl}}} \sqrt{\frac{\gamma}{r}} M \left(1 + 0.5 \gamma_1 M^2\right)^{-\frac{\gamma+1}{2(\gamma-1)}} Y. \quad (A6)$$
For the series of flows we consider, with $M = 1$ at the throat and a fixed stagnation state at the inlet, a consequence of Equation (A6) is that the only way to increase the mass flow is to increase the throat section. This was achieved by the performed optimization, which increased the nozzle throat and the mass flow and, as a consequence, the thrust (more precisely, the mass flow was increased by 4.19% with respect to its nominal value at the first optimization step and then by an additional 0.78% at the second step).

## References

1. Ferri, A. *Application of the Method of Characteristics to Supersonic Rotational Flow*, TR 841; NASA: Washington, DC, USA, 1946.

2. Shapiro, A. *The Dynamics and Thermodynamics of Compressible Fluid Flow*; The Ronald Press Company: New York, NY, USA, 1954.

3. Bonnet, A.; Luneau, J. *Aérodynamique. Théories de la Dynamique des Fluides*; Cepadues: Toulouse, France, 1989.

4. Délery, J. *Traité D'aérodynamique Compressible*; Collection Mécanique des Fluides, Lavoisier—Hermés, Science: Cachan, France, 2008; Volume 3.

5. Jameson, A. Aerodynamic design via control theory. *J. Sci. Comput.* **1988**, *3*, 233–260. [CrossRef]

6. Brezillon, J.; Gauger, N. 2D and 3D Aerodynamic Shape Optimization Using the Adjoint Approach. *Aerosp. Sci. Technol.* **2004**, *8*, 715–727. [CrossRef]

7. Brezillon, T.; Dwight, R. Discrete adjoint of the Navier-Sotkes equations for aerodynamic shape optimization. In *Proceedings of the EUROGEN 2005*, Munich, Germany, 12–14 September 2005.

8. Peter, J.; Dwight, R. Numerical sensitivity analysis for aerodynamic optimization: A survey of approaches. *Comput. Fluids* **2010**, *39*, 373–391. [CrossRef]

9. Schmidt, S.; Ilic, C.; Schultz, V.; Gauger, N. Three dimensional large scale aerodynamic shape optimization based on shape calculus. *AIAA J.* **2013**, *51*, 2615–2627. [CrossRef]

10. Venditti, D.; Darmofal, D. Grid adaptation for functional outputs: Application to two-dimensional inviscid flows. *J. Comput. Phys.* **2002**, *176*, 40–69. [CrossRef]

11. Dwight, R. Heuristic a *posteriori* estimation of error due to dissipation in finite volume schemes and application to mesh adaptation. *J. Comput. Phys.* **2008**, *227*, 2845–2863. [CrossRef]

12. Todarello, G.; Vonck, F.; Bourasseau, S.; Peter, J.; Désidéri, J.A. Finite-volume goal-oriented mesh-adaptation using functional derivative with respect to nodal coordinates. *J. Comput. Phys.* **2016**, *313*, 799–819. [CrossRef]

13. Giles, M.; Duta, M.; Müller, J.D. Adjoint code developments using the exact discrete approach. In *Proceedings of the 15th AIAA Computational Fluid Dynamics Conference, Anaheim, CA, USA, 11–14 June 2001; AIAA Paper Series, Paper 2001-2596*.

14. Giles, M.; Duta, M.; Müller, J.D.; Pierce, N. Algorithm developments for discrete adjoint methods. *AIAA J.* **2003**, *41*, 198–205. [CrossRef]

15. Economon, T.; Alonso, J.; Albring, T.; Gauger, N. Adjoint Formulation Investigations of Benchmark Aerodynamic Design Cases in SU2. In *Proceedings of the 35th AIAA Applied Aerodynamics Conference, Denver, CO, USA, 5–9 June 2017; AIAA Paper Series, Paper 2017-4363*.

16. Sagebaum, M.; Albring, T.; Gauger, N. High-Performance Derivative Computations using CoDiPack. *ACM Trans. Math. Softw.* **2019**, *45*, 1–26. [CrossRef]

17. Pinel, X.; Montagnac, M. Block Krylov methods to solve adjoint problems in aerodynamic design optimization. *AIAA J.* **2013**, *51*, 2183–2189. [CrossRef]

18. Xu, S.; Radford, D.; Meyer, M.; Müller, J.D. Stabilisation of discrete steady adjoint solvers. *J. Comput. Phys.* **2015**, *299*, 175–195. [CrossRef]

19. Lions, J.L. *Contrôle Optimal de Systèmes Gouvernés par des Équations aux Derivées Partielles*; Etudes Mathématiques; Dunod, Gauthier-Villars: Paris, France, 1968.

20. Sartor, F.; Mettot, C.; Sipp, D. Stability, Receptivity, and Sensitivity Analyses of Buffeting Transonic Flow over a Profile. *AIAA J.* **2015**, *53*, 1980–1993. [CrossRef]

21. Laurenceau, J. Surfaces de réponses par Krigeage pour l'optimisation de formes. Ph.D. Thesis, Institut National Polytechnique de Toulouse, Toulouse, France, 2008.

22. de Baar, J.; Scholcz, T.; Dwight, R. Technical Note: Exploiting Adjoint Derivatives in High-Dimensional Metamodels. *AIAA J.* **2015**, *53*, 1391–1395. [CrossRef]

23. Sipp, D.; Marquet, O.; Meliga, P.; Barbagallo, A. Dynamics and Control of Global Instabilities in Open-Flows: A Linearized Approach. *Appl. Mech. Rev.* **2010**, *63*, 030801. [CrossRef]

24. Luchini, P.; Bottaro, A. Adjoint Equations in Stability Analysis. *Annu. Rev. Fluid Mech.* **2014**, *46*, 493–517. [CrossRef]

25. Parish, E.; Duraisamy, K. A paradigm for data-driven predictive modeling using field inversion and machine learning. *J. Comput. Phys.* **2016**, *305*, 758–774. [CrossRef]

26. Talagrand, O.; Courtier, P. Variational Assimilation of Meteorological Observations With the Adjoint Vorticity Equation. I: Theory. *Q. J. R. Meteorol. Soc.* **1987**, *113*, 1311–1328. [CrossRef]

27. Hirsch, C. *Numerical Computation of Internal and External Flows: The Fundamentals of Computational Fluid Dynamics*, 2nd ed.; Butterworth—Heineman; Elsevier: Amsterdam, The Netherlands, 2007.
28. Anderson, W.; Venkatakrishnan, V. Adjoint design optimization on unstructured grid using a continuous adjoint formulation. *Comput. Fluids* **1999**, *28*, 443–480. [CrossRef]

29. Giles, M.; Pierce, N. Analytic adjoint solutions for the quasi-one-dimensional Euler equations. *J. Fluid Mech.* **2001**, *426*, 327–345. [CrossRef]

30. Peter, J.; Désidéri, J.A. Ordinary differential equations for the adjoint Euler equations. *Phys. Fluids* **2022**, *34*, 086113. [CrossRef]

31. Meyer, R.; Goldstein, S. The method of characteristics for problems of compressible flow involving two independent variables: Part I. The general theory. *Q. J. Mech. Appl. Math.* **1948**, *1*, 196–219. [CrossRef]

32. Anderson, J. *Modern Compressible Flow*, 3rd ed.; McGraw-Hill Series in Aeronautical and Aerospace Engineering; McGraw-Hill: New York, NY, USA, 2003.

33. Oberkampf, W.; Trucano, T. Verification and validation in computational fluid dynamics. *Prog. Aerosp. Sci.* **2002**, *38*, 209–272. [CrossRef]

34. Liepmann, H.; Roshko, A. *Elements of Gasdynamics*; Dover Publications: Mineola, NY, USA, 1956.

35. Grossman, B. *Fundamental Concepts of Real Gasdynamics*; Virginia Polytechnic Institute and State University: Blacksburg, VA, USA, 2000.

36. Mourouzidis, C. Preliminary Design of Next Generation Mach 1.6 Supersonic Business Jets to Investigate Landing and Take-Off (LTO) Noise and Emissions. In Proceedings of the 12th EASN conference on Innovation and Space for opening New Horizons, Barcelona, Spain, 18–21 October 2022.

37. Cambier, L.; Heib, S.; Plot, S. The elsA CFD software: Input from research and feedback from industry. *Mech. Ind.* **2013**, *14*, 159–174. [CrossRef]

38. Jameson, A.; Schmidt, W.; Turkel, E. Numerical Solutions of the Euler Equations by Finite Volume Methods Using Runge-Kutta Time-Stepping Schemes. In Proceedings of the 14th Fluid and Plasma Dynamics Conference, Palo Alto, CA, USA, 23–25 June 1981; AIAA Paper Series, Paper 1981-1259.

39. Peter, J.; Renac, F.; Dumont, A.; Méheut, M. Discrete adjoint method for shape optimization and mesh adaptation in the elsA code. *Status and Challenges. In Proceedings of the 50th 3AF Symposium on Applied Aerodynamics, Toulouse, France, 29 March–1 April 2015.*

40. Peter, J.; Renac, F.; Labbé, C. Analysis of finite-volume discrete adjoint fields for two-dimensional compressible Euler flows. *J. Comput. Phys.* **2022**, *449*, 110811. [CrossRef]

41. Giles, M.; Pierce, N. Adjoint equations in CFD: Duality, boundary conditions and solution behaviour. In Proceedings of the 13th Computational Fluid Dynamics Conference, Snowmass Village, CO, USA, 29 June–2 July 1997; AIAA Paper Series, Paper 97-1850.

42. Lozano, C. On mesh sensitivities and boundary formulas for discrete adjoint-based gradients in inviscid aerodynamic shape optimization. *J. Comput. Phys.* **2017**, *346*, 403–436. [CrossRef]

43. Lozano, C. Singular and discontinuous solutions of the adjoint Euler equations. *AIAA J.* **2018**, *56*, 4437–4451. [CrossRef]

44. Lozano, C. Watch your adjoints! Lack of mesh convergence in inviscid adjoint solutions. *AIAA J.* **2019**, *57*, 3991–4006. [CrossRef]

45. Peter, J. *Contributions to Discrete Adjoint Method in Aerodynamics for Shape Optimization and Goal-Oriented Mesh Adaptation*; University of Nantes. Mémoire pour Habilitation à Diriger des Recherches: Nantes, France, 2020.

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.