Analysis of finite-volume discrete adjoint fields for two-dimensional compressible Euler flows
AAnalysis of finite-volume discrete adjoint fields for two-dimensionalcompressible Euler flows
Jacques Peter a, ∗ , Florent Renac a , Clément Labbé a a DAAA, ONERA, Université Paris Saclay, F-92322 Châtillon, France
Abstract
This work deals with a number of questions relative to the discrete and continuous adjoint fields associated to Eulerequations and classical aerodynamic functions. The consistency of the discrete adjoint equations with the correspond-ing continuous adjoint partial differential equation is one of them. It is has been established or at least discussed onlyfor a handful of numerical schemes and a contribution of this article is to give the adjoint consistency conditions forthe 2D Jameson-Schmidt-Turkel scheme in cell-centred finite-volume formulation. The consistency issue is also stud-ied here from a new heuristic point of view by discretizing the continuous adjoint equation for the discrete flow andadjoint fields. Both points of view prove to provide useful information. Besides, it has been often noted that discreteor continuous inviscid lift and drag adjoint exhibit numerical divergence close to the wall and stagnation streamlinefor a wide range of subsonic and transonic flow conditions. This is analyzed here using the physical source term per-turbation method introduced in reference [1]. With this point of view, the fourth physical source term of [1] appearsto be the only one responsible for this behavior. It is also demonstrated that the numerical divergence of the adjointvariables corresponds to the response of the flow to the convected increment of stagnation pressure and diminution ofentropy created at the source and the resulting change in lift and drag.
Keywords: discrete adjoint method, continuous adjoint method, dual consistency, adjoint Rankine-Hugoniot,compressible Euler equations
1. Introduction
Discrete and continuous adjoint methods are well established methods to efficiently calculate derivatives of aero-dynamic functions with respect to numerous design parameters. Adjoint-based derivatives and adjoint-fields arecommonly used for local shape optimization [2, 3, 4, 5, 6], goal-oriented mesh adaptation [7, 8, 9, 10, 11, 12, 13, 14],and also other fields like flow control [15, 16], meta-modelling [17], receptivity-sensitivity-stability analyses [18] anddata assimilation [19]. These methods are often used for the linear analysis of non-linear conservation laws where theadjoint is defined as the dual to the linearized equations around a given solution to the direct non-linear problem. Thecontinuous adjoint method refers to the discretization of the adjoint equations associated to the formal direct problem,while the discrete adjoint method refers to the adjoint equations to the discrete direct problem.Concerning the discrete adjoint method, strong efforts have been devoted to implement (either by hand or usingautomatic differentiation) the exact linearization of complex schemes for complex models [20, 21] and also to enhancethe robustness of the linear solver [22, 23] whereas fundamental questions relative to discrete adjoint fields are stillopen even for two-dimensional (2D) compressible Euler flows. Numerical simulations including viscous laminar andturbulent effects are today common practices in CFD. It is nevertheless desirable to gain understanding in the discreteadjoint method addressing unresolved issues in the framework of 2D Euler equations.First, are the discrete adjoint fields consistent with the continuous equation at the limit of fine meshes? This ques-tion goes back to a discussion by Giles et al. in [24] where the authors discuss the discrete adjoint counterpart of a ∗ Corresponding author. Tel.: +33 1 46 73 41 84.
Email addresses: [email protected] (Jacques Peter), [email protected] (Florent Renac)
Preprint submitted to Elsevier September 16, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] S e p trong boundary condition at a wall and derive a normal projection of the momentum adjoint variable that does notcorrespond to the normal momentum component of the continuous adjoint solution at the boundary. This of courseraises the question which schemes are adjoint-consistent and which are not, and what are the consequences of the lackof dual-consistency. These questions have been studied in a handful of articles where the authors have first to preciselydescribe the numerical method they consider, including the possibly specific discretization close to the boundary, andalso the way the function of interest is discretized. Lu and Darmofal [25] considered a discontinuous Galerkin (DG)scheme and then demonstrated how to formulate the primal boundary condition and the output functional to obtaindual-consistency. The difference between the consistent and non-consistent formulations is noticed in both the appear-ance of the adjoint contours close to the function support and the convergence rate of the output of interest w.r.t. thecharacteristic mesh size. Hartmann [26] also discussed dual consistency for DG schemes, including interior penaltyDG and also dealt with Navier-Stokes equations and test-cases. The adjoint consistency of a high-order correctionprocedure via reconstruction formulation for hyperbolic conservation laws was investigated in [27]. As concerningfinite-differences and finite-volume (FV), the questions of adjoint-consistency was discussed by Duivesteijn et al.and Lozano for quasi 1D Euler flows [28, 29] and also for an unsteady conservation law by Liu and Sandu [30]. Inthis later article, the authors considered first-order and second-order upwind schemes and derived the circumstancesin which the discrete adjoint is punctually not consistent (change of convection direction, inadequate discretizationof boundary condition in particular). Hicken and Zingg have demonstrated the influence of adjoint consistency for aclass of summation-by-parts difference schemes, illustrating the benefit of consistency in terms of regularity of adjointcontours close to boundaries, functional mesh convergence and dual-weighted corrected functional mesh convergence[31]. Recently Stück carried out a corresponding effort for cell-vertex FV discretization (using median dual grids),deriving the condition on the output of interest and scheme residual for dual-consistency [32, 33]. In this article, theconditions for dual consistency of the classical Jameson-Schmidt-Turkel (JST) scheme in 2D cell-centered FV arederived. Besides, the dual consistency is also studied here from a heuristic point of view, discretizing the continuousadjoint equation for the discrete adjoint fields which adds valuable information to the theoretical study where the flowor the gradient of the adjoint is discontinuous and where numerical divergence of the adjoint is observed.Then, in the case of hyperbolic equations, the validity of the adjoint linearization around discontinuities in the directsolution has to be carefully handled because it results in linear equations with discontinuous coefficients. The analysismust include the linearization of the jump relations at the discontinuity [34] which leads to a so-called interior bound-ary condition for the adjoint variables [35]. The adjoint relations to the Rankine-Hugoniot (RH) relations have beenderived in [36] for the quasi-one-dimensional and in [37, 38] for the 2D compressible Euler equations and provedcontinuity of the adjoint variables across the shock. In contrast, their derivatives may be discontinuous [36, 1, 37].Lozano [38] also briefly discussed the consistency of the continuous adjoint fields but on rather coarse meshes. Theconclusions of his work leave the question open whether these relations are actually satisfied at the limit of fine dis-crete adjoint fields.Finally, it has been well-documented [1, 39, 38, 40] that, at certain flow conditions, the inviscid lift and drag adjointfields may exhibit increasing values as the mesh is refined close to the stagnation streamline, the wall and, for (atleast locally) supersonic flows, the Mach lines impacting a shock foot. (An interested reader may also consider thesame problem for a 2D potential flow [41, 40].) In the case of the stagnation line, Giles and Pierce derived a 1 / √ d law for lift and drag adjoint [1] ( d being the distance to the line), but they used the potential flow theory and this1 / √ d behaviour is not always well-observed for compressible flows [39, 38]. In an effort to analyse the zones ofsingular behavior, Todarello et al. [42] used the characterization of the discrete adjoint and looked how the flow isperturbed by a residual perturbation located in one of the these specific areas. This method is used again here but withperturbations of the residual corresponding to physical source terms [1] which allows fluid mechanics analysis. In thisframework, it appears that the numerical divergence of the adjoint variables is due to the influence of the fourth of thesource term proposed in [1], the corresponding convected increment of stagnation pressure and diminution of entropy(created at the source location) and the resulting change in lift and drag. Finally, another aspect of lift and drag 2Dinviscid adjoint field is studied in details: in the case of a supersonic flow with detached shock wave, the structure ofthe adjoint field is theoreticaly derived in the constant flow area and the consistency of fine grid numerical solutionswith this analytical model is checked.The paper is organized as follows. Sections 2 and 3 are devoted to reminders about continuous and discrete adjointsolutions for Euler flows. The adjoint consistency of the cell-centered JST scheme with structured 2D meshes [43] isdiscussed in § 4. Section 5 presents the adjoint RH equations for Euler flows. The theoretical results about adjoint2onsistency and adjoint-gradient discontinuity at shock waves are then assessed in § 6 using fine grid lift- and drag-adjoint fields for inviscid flows around the NACA0012 airfoil in supersonic (§ 6.1), transonic (§ 6.2) and subsonic(§ 6.3) regimes. In these sections, the question of the asymptotic behavior of adjoint at the wall/stagnation line andmore generally high gradient/high-values in the adjoint fields is also discussed using Giles and Pierce physical sourceterm approach derived in § 3.3. Finally, concluding remarks about this work are given in § 7.
2. Continuous adjoint equations for 2D Euler flows
The continuous adjoint equations for compressible flows were first derived by Jameson [4] in the case of a 2DEuler flow about a profile. He considered a body fitted structured grid that was mapped to a conformal rectangle andused the Euler equations in the resulting ( ξ , η ) -coordinates. A parametrization of the mapping then allowed to varythe airfoil shape in the physical space (without altering the domain of variation of the transformed coordinates) and todefine a gradient calculation problem for functional outputs. The adjoint were also derived in [44, 36, 38, 37] using aframework including the linearization of the RH relations.As, in this formulation, the system of transformed coordinates is attached to a structured mesh, the aforementionedequations could not be used for unstructured CFD for which a formulation in physical coordinates was necessary. Thecorresponding system of equations was first derived by Anderson and Venkatakrishnan in [45] (and one year laterby Hiernaux and Essers [46, 47]). A slightly simplified presentation of the theoretical part of reference [45] is givenbelow.The quantity of interest is assumed to be the projection of the force applied by the fluid onto the solid, projectedin the direction d = ( − cos α , sin α ) T for the lift, or d = ( sin α , cos α ) T for drag with α the angle of attack: J = (cid:90) Γ w p ( n · d ) ds , (1)where Γ w is the boundary of the solid body, n = ( n x , n y ) T is the local normal (external for fluid and internal for thesolid) and p is the static pressure.The 2D compressible Euler equations in conservative form read ∂ F x ( W ) ∂ x + ∂ F y ( W ) ∂ y = Ω , (2) W being the vector of conservative variables and F x and F y the physical fluxes, W = ρρ u x ρ u y ρ E , F x ( W ) = ρ u x ρ u x + p ρ u x u y ρ u x H , F y ( W ) = ρ u y ρ u x u y ρ u y + p ρ u y H , with H = E + p / ρ the total specific enthalpy. We consider an ideal gas law for the static pressure p = ( γ − ) ρ e = ( γ − )( ρ E − ρ (cid:107) U (cid:107) ) , with γ the ratio of specific heats and U = ( u x , u y ) T the velocity vector. Let δ W be a perturbationin the steady state flow that is caused be an infinitesimal perturbation of airfoil shape or flow conditions. As W and W + δ W are solutions of the steady Euler equations for initial and perturbed problem, by difference we get ∂ ( A δ W ) ∂ x + ∂ ( B δ W ) ∂ y = Ω , with A and B the Jacobians of the fluxes F x and F y , respectively. The perturbation in J can be augmented by the dotproduct of above equation with an arbitrary co-state field ψδ J = (cid:90) Γ w δ p ( n · d ) ds + (cid:90) Γ w p ( δ ( n ) · d ) ds + (cid:90) Γ w p ( n · d ) δ ( ds ) + (cid:90) Ω ψ T (cid:18) ∂ ( A δ W ) ∂ x + ∂ ( B δ W ) ∂ y (cid:19) dv . Assuming smooth direct and adjoint solutions, the last term can be transformed by integration by parts into − (cid:90) Ω (cid:18) ∂ ψ T ∂ x A + ∂ ψ T ∂ y B (cid:19) δ W dv + (cid:90) Γ w ψ T ( n x A + n y B ) δ W ds + (cid:90) Γ ∞ ψ T ( n x A + n y B ) δ W ds , Γ ∞ the farfield boundary. The perturbation in the objective may then be rewritten as δ J = (cid:90) Γ w δ p ( n · d ) ds + (cid:90) Γ w p ( δ ( n ) · d ) ds + (cid:90) Γ w p ( n · d ) δ ( ds ) − (cid:90) Ω (cid:18) ∂ ψ T ∂ x A + ∂ ψ T ∂ x B (cid:19) δ W dv + (cid:90) Γ w ψ T ( n x A + n y B ) δ W ds + (cid:90) Γ ∞ ψ T ( n x A + n y B ) δ W ds (3)The adjoint method removes the dependency in the flow perturbation δ W for the calculation of the variation of J .This directly yields the adjoint equation in the fluid domain − A T ∂ ψ∂ x − B T ∂ ψ∂ y = , in Ω . (4)Note that the last equation in (4) reads U · ∇ ψ − HU · ∇ ψ = . (5)Besides the wall boundary condition on Γ w links δ n and δ W : δ ( U · n ) = δ ( U ) · n + U · δ ( n ) = δ W satisfy δ W n x + δ W n y + W δ n x + W δ n y = . (6)The explicit calculation of ψ T ( n x A + n y B ) δ W at the wall yields ψ T ( n x A + n y B ) δ W = ( n x ψ + n y ψ ) (cid:18) ( γ − ) (cid:107) U (cid:107) , ( − γ ) u x , ( − γ ) u y , ( γ − ) (cid:19) δ W δ W δ W δ W + ( ψ + u ψ + v ψ + H )( n x δ W + n y δ W ) where the first term of the right-hand-side accounts for the differentiated pressure term of Euler flux and the secondterm for the differentiated convective part of Euler flux. The last term is reformulated using the linearized boundarycondition (6) then the sum of the terms involving the flow variation δ W in the integrals over Γ w is set to zero to definethe adjoint wall boundary condition ( n · d ) d pdW δ W + ( n x ψ + n y ψ ) (cid:18) ( γ − ) (cid:107) U (cid:107) , ( − γ ) u x , ( − γ ) u y , ( γ − ) (cid:19) δ W δ W δ W δ W = . (7)The last term of (7) corresponds to the derivative of the static pressure with respect to W , so that the final andclassical form of the boundary condition is simply n · d + ψ n x + ψ n y = . (8)In the farfield, no variation of the boundary needs to be considered. The Jacobian in direction n can be rewrittenby using a locally one dimensional characteristic decomposition to yield (cid:90) Γ ∞ ψ T ( n x A + n y B ) δ W ds = (cid:90) Γ ∞ ( ψ T P − DP δ W ) ds . (9)It is assumed that P δ W (cid:39) δ ( PW ) . The variation in these characteristic variables is zero for the componentscorresponding to negative eigenvalues of the Jacobian, n x A + n y B , in the classical 1D approximate linearization at theboundary (information coming from outside of the domain and fixed characteristic value). The farfield adjoint BCsimply imposes that the other component of ψ T P − are zero so that δ W T ( n x A T + n y B T ) ψ vanishes for all δ W .4 . Discrete adjoint equations The finite volume scheme of interest defines the steady-state discrete flow vector W (of size n W ) as the solution ofa set of n W non-linear equations involving W and the vector of mesh coordinates X : R ( W , X ) = . Let us assume that the mesh is a regular function of a vector of design parameters β (of size n β ), then the implicitfunction theorem allows to define W as a function of X and so of β [48]. Discrete gradient calculation consists incomputing the derivatives of n f functions J k ( β ) = J k ( W ( β ) , X ( β )) , ≤ k ≤ n f , with respect to the n β design parameters. In external aerodynamic applications, n β is usually much larger than n f andthe most efficient way to proceed is to use the discrete adjoint method which requires to solve n f linear systems (cid:18) ∂ R ∂ W (cid:19) T Λ k = − (cid:18) ∂ J k ∂ W (cid:19) T , ≤ k ≤ n f , (10)then to calculate d J k d β = ∂ J k ∂ X d X d β + Λ Tk (cid:18) ∂ R ∂ X d X d β (cid:19) . The dominant cost is the inversion of the n f linear systems of size n W whereas all other classical methods solve n β linear (or possibly non-linear) systems of size n W . The adjoint vector Λ associated to one function J is usually identified to the sensitivity of J to a perturbationof the residual R followed by re-convergence (see [12] and references therein). Following [1], we consider here aperturbation δ R added in the right-hand-side of the discrete flow equations and denote by W + δ W the convergedsolution corresponding to the perturbed equation: R ( W + δ W , X ) = δ R , or at first order − δ R + ∂ R ∂ W δ W = . Since J ( W + δ W , X ) (cid:39) J ( W , X ) + ∂ J ∂ W δ W , the first order change in the function of interest J due to change in flow δ W is δ J = ∂ J ∂ W (cid:16) ∂ R ∂ W (cid:17) − δ R . Involving the discrete adjoint vector Λ , yields δ J = − Λ T δ R . (11)If only the a -th component of R at cell m has been arbitrarily altered by a small number δ R am , then (11) yields Λ am = − δ J / δ R am , (12)which defines the a -th component of Λ at cell m as the limit ratio of the change in J divided by the infinitesimal changein the residual R at the corresponding cell and component which caused the change in flow and function value.5 .3. Physical characterization for 2D Euler flows We here recall the physical interpretation of the adjoint solution by Giles and Pierce [1]. More precisely, thediscussion in reference [1] is based on the continuous adjoint equation and its discrete counterpart is presented here.Four physical source terms δ R are defined at each individual cell: 1) local mass source at fixed stagnation pressureand total enthalpy; 2) local normal force; 3) local change in total enthalpy at fixed static and total pressure; 4) localchange in total pressure at fixed total enthalpy and static pressure: δ R m = ε u x u y H , δ R m = ε − ρ u y ρ u x , δ R m = ε − H , δ R m = ε (cid:16) γ − γ + γ M (cid:17) u x (cid:16) γ − γ + γ M (cid:17) u y (cid:16) γ − γ + γ M (cid:17) H (cid:16) γ − γ + γ M (cid:17) , < ε (cid:28) , (13) where p , H and M denote the stagnation pressure, total enthalpy and Mach number, respectively.At first order, the corresponding functions variations δ J m , δ J m , δ J m , δ J m due to the source terms in cell m , (13 ),read ( δ J m , δ J m , δ J m , δ J m ) = − ( Λ m , Λ m , Λ m , Λ m ) × δ R , m δ R , m δ R , m δ R , m δ R , m δ R , m δ R , m δ R , m δ R , m δ R , m δ R , m δ R , m δ R , m δ R , m δ R , m δ R , m , (14)where δ J dm corresponds to the change in J due to d -th change of R in cell m . Reconverging flow W with a globalchange δ W d leads to the change δ J dm of the function of interest. As the four changes in R are linearly independent,Giles and Pierce define the adjoint vector at cell m as the solution to (14) so, using the expression of the perturbationterms (13), we obtain ( Λ m , Λ m , Λ m , Λ m ) = − ( δ J m , δ J m , δ J m , δ J m ) × ε − − ( γ − ) M − ( +( γ − ) M ) u x (cid:107) U (cid:107) − ( +( γ − ) M ) u y (cid:107) U (cid:107) + ( γ − ) M H − u y ρ (cid:107) U (cid:107) u x ρ (cid:107) U (cid:107) − H − γ M γ u x c γ u y c − γ M H . (15)The residual perturbation δ R has thus been physically defined which makes the local discrete adjoint vectorphysically defined. For all systems of equations for which a similar demonstration can be done, the adjoint vector getsintrinsic and we can expect similar solutions from different discretizations and also mesh convergence. Let us finallynote that not only the Λ components is of physical significance, but also the ( Λ · δ R d ) / ε terms.Finally, Giles and Pierce [1] have noted for continuous adjoint that the third perturbation does not alter the pressurefield for inviscid flows and hence lets the drag and lift unchanged. From the expression of this perturbation, it is thenstraighforward to prove that lift and drag continuous adjoint fields satisfy all over the fluid domain the followingequation: ψ = H ψ . (16)This property may be compared to (5) and is seen to be well satisfied in the numerical adjoint fields (see § 6.2.2).The fourth perturbation δ R also preserves the static pressure at the source term location but the authors of refer-ence [1] prove that, contrary to δ R , for arbitrary flow conditions, it does preserve neither the static pressure fieldnor the streamtube structure of the base flows (in particular, the mass flux varies downwind the source term in thecorresponding streamtube of the original flow). 6 . Adjoint-consistency of the cell-centered JST scheme The question of dual consistency of the 2D cell-centered FV JST scheme is here analyzed by using the classicalmethod of equivalent differential equation which consists in introducing Taylor series developments for all ajointvalues to analyze the truncation error of the scheme [49, ch. 9]. We stress that this approach remains valid in the caseof discontinuous solutions [50]. This analysis will be completed in § 6.1 to § 6.3 by the examination of the residualsof the continuous adjoint equation calculated for the discrete adjoint field. In the following we assume that the JSTscheme is differentiable w.r.t. the degrees of freedom of the discrete direct problem which is a necessary condition foradjoint consistency. Situations of non-differentiability are discussed at the end of § 4.1.Lozano studied the adjoint consistency of the discrete adjoint of the JST scheme for both cell-centered and cell-vertex discretizations [29] of the quasi-one-dimensional compressible Euler equations. Boundary conditions usingartifical dissipation are applied which lead to unconsistent discrete adjoint scheme at the boundaries. Such boundaryconditions are also not suitable for the simulation of internal flows that require conservativity. In this work wecircumvent these difficulties by removing the dissipation at the boundary faces, while we modify the linearization atthe penultimate cell to wall boundaries.Jameson presented the continuous adjoint equations of 2D and 3D Euler equations in transformed coordinates [4].Later on, Giles and Pierce derived the corresponding direct differentiation equation in 2D [51]. In this framework, thephysical space (typically about an airfoil) is seen a the image of the unit domain [ , ] × [ , ] with coordinates ( ξ , η ) .The transformed 2D Euler equations in the ( ξ , η ) -coordinates are well-known and, for functions of the form (1), theadjoint equations read − ∂ Λ T ∂ ξ (cid:18) ∂ y ∂ η A ( w ) − ∂ x ∂ η B ( w ) (cid:19) − ∂ Λ T ∂ η (cid:18) − A ( w ) ∂ y ∂ ξ + B ( w ) ∂ x ∂ ξ (cid:19) = , (17) A ( w ) and B ( w ) being the Jacobians of Euler flux in the physical x and y directions. The wall being defined by the ξ = d = ( d x , d y ) , thenthe corresponding boundary condition − Λ T (cid:18) A ( w ) ∂ y ∂ η − B ( w ) ∂ x ∂ η (cid:19) − ( ∂ y ∂ η d x − ∂ x ∂ η d y ) ∂ p ∂ W = ξ = Λ ∂ y ∂ η − Λ ∂ x ∂ η + ( ∂ y ∂ η d x − ∂ x ∂ η d y ) = . (18)The adjoint consistency of the FV scheme of interest is discussed hereafter referring to the analytic equations (17)and (18).Using a structured mesh, the JST scheme in the i -direction reads (the j subscripts have been dropped for the sakeof readability): F JSTi + / = ( F ( W i ) + F ( W i + )) · S i + / − k ν i + / κ i + / ( W i + − W i )+ k i + / κ i + / ( W i + − W i + + W i − W i − ) (19) ν i = | p i + − p i + p i − | ( p i + + p i + p i − ) ν i + / = max ( ν i , ν i + ) (20) k i + / = max ( , k − ν i + / k ) κ i + / = | U i + / . S i + / | + c i + / (cid:107) S i + / (cid:107) , (21)the fluid velocity and sound velocity at interface i + / F JST / = ( F ( W ) + F ( W )) · S / − k ν / κ / ( W − W ) + k / κ / ( W − W + W ) (22) F JST / = F ( W b ) · S / (23)7here p = p b − p is used in the evaluation of ν , and W b is a boundary state satisfying all Dirichlet-like boundaryconditions.Of course, at specific faces where either U · S =
0, or k − ν k =
0, or ν i = ν i + , the flux is not differentiable. Thisissue is discussed in [52] for a simulation about a symmetric airfoil, with a symmetric mesh and zero angle of attackleading to two complete lines of faces where | U · S | is extremely close to machine zero. The small irregularities ofadjoint quantities that appear close to these lines are cured by regularization but no such problem is observed for morecomplex cases and, most generally, this question is discarded [53] considering that the occurrence of such equalitiesis marginal in practice.Finally, we assume that the increments ∆ ξ and ∆ η between successive mesh points, have the same order ofmagnitude. Due to the four-point single-directional stencil, the discrete adjoint equation at a current cell ( i , j ) of the structuredmesh reads ( Λ T ( i − , j ) − Λ T ( i − , j ) ) ∂ F i − / , j ∂ W ( i , j ) + ( Λ T ( i − , j ) − Λ T ( i , j ) ) ∂ F i − / , j ∂ W ( i , j ) (24) +( Λ T ( i , j ) − Λ T ( i + , j ) ) ∂ F i + / , j ∂ W ( i , j ) + ( Λ T ( i + , j ) − Λ T ( i + , j ) ) ∂ F i + / , j ∂ W ( i , j ) (25) +( Λ T ( i , j − ) − Λ T ( i , j − ) ) ∂ F i , j − / ∂ W ( i , j ) + ( Λ T ( i , j − ) − Λ T ( i , j ) ) ∂ F i , j − / ∂ W ( i , j ) (26) +( Λ T ( i , j ) − Λ T ( i , j + ) ) ∂ F i , j + / ∂ W ( i , j ) + ( Λ T ( i , j + ) − Λ T ( i , j + ) ) ∂ F i , j + / ∂ W ( i , j ) = . (27)Only fluxes in the i -direction appear in lines (24)-(25), while only fluxes in the j -direction appear in lines (26)-(27). It is possible to analyze separately the two groups of corresponding terms. This is done below in the i -direction(again dropping the j subscripts). The contribution of the centered flux in the derivation of (24)-(25) is (cid:16) ( Λ Ti − − Λ Ti ) ∂ F ∂ W i S i − / + ( Λ Ti − Λ Ti + ) ∂ F ∂ W i S i + / (cid:17) . (28)The i and j subscripts and ξ and η coordinates are linked by simple affine transformations ξ ( i ) = i − i max − = ( i − ) ∆ ξ η ( j ) = j − j max − = ( j − ) ∆ η and the expression of the surface vectors is easily derived from those of the coordinates S i + / , j = (cid:18) y ( ξ ( i + / ) , η ( j + / )) − y ( ξ ( i + / ) , η ( j − / )) − x ( ξ ( i + / ) , η ( j + / )) + x ( ξ ( i + / ) , η ( j − / )) (cid:19) . The two terms in (28) are hence consistent with − ∂ Λ T ∂ ξ (cid:18) ∂ y ∂ η A − ∂ x ∂ η B (cid:19) ∆ ξ ∆ η respectively at points ( i − / , j ) and ( i + / , j ) , with ∆ ξ ∆ η denoting the volume of a cell, and the sum of bothcontributions is a consistent approximation of − ∂ Λ T ∂ ξ (cid:18) ∂ y ∂ η A − ∂ x ∂ ξ B (cid:19) ∆ ξ ∆ η at point ( i , j ) . Moreover it is second order accurate as it uses symmetric means and differences w.r.t. ( ξ ( i ) , η ( j )) .Carrying out the corresponding expansions for the linearized fluxes in j -direction, the second term of equation (17)8imes ∆ ξ ∆ η is recovered. The first part of the discrete adjoint equation, obtained by isolating the derivation of thecentered flux, is hence a consistent discretization of equation (17) up to a multiplicative ∆ ξ ∆ η factor.After this first step, the discussion of the interior-cell adjoint consistency of JST scheme comes back to the exam-ination of all the other terms of the equation (24)-(27). and the question whether other second order terms in ∆ ξ and ∆ η appear (no consistency) or only higher order terms appear (consistency). The terms of (24)-(25) that involve thederivation of the sensor, read k (cid:18) ( Λ Ti − − Λ Ti − ) ∂ ν i − / ∂ W ( i , j ) κ i − / (cid:16) − ( W i − − W i − ) − [ k i − / > ] ( W i − W i − + W i − − W i − ) (cid:17) +( Λ Ti − − Λ Ti ) ∂ ν i − / ∂ W ( i , j ) κ i − / (cid:16) − ( W i − − W i ) − [ k i − / > ] ( W i + − W i + W i − − W i − ) (cid:17) +( Λ Ti − Λ Ti + ) ∂ ν i + / ∂ W ( i , j ) κ i + / (cid:16) − ( W i − W i + ) − [ k i + / > ] ( W i + − W i + + W i − W i − ) (cid:17) +( Λ Ti + − Λ Ti + ) ∂ ν i + / ∂ W ( i , j ) κ i + / (cid:16) − ( W i + − W i + ) − [ k i + / > ] ( W i + − W i + + W i + − W i ) (cid:17)(cid:19) , (29)where 1 [ k i + / > ] stands for 1 if k i + / is strictly positive and 0 if not. As the spectral radius κ is O ( ∆ η ) , all the termsin (29) are at least O ( ∆ ξ ∆ η ) . The terms of (24)-(25) that involve the derivation of the spectral radius are ( Λ Ti − − Λ Ti ) ∂ κ i − / ∂ W ( i , j ) (cid:16) − k ν i − / ( W i − − W i ) + k i − / ( W i + − W i + W i − − W i − ) (cid:17) ( Λ Ti − Λ Ti + ) ∂ κ i + / ∂ W ( i , j ) (cid:16) − k ν i + / ( W i − W i + ) + k i + / ( W i + − W i + W i − W i − ) (cid:17) . (30)As κ i − / and its derivatives w.r.t. W are O ( ∆ η ) and as ν i − / is O ( ∆ ξ ) these terms are O ( ∆ ξ ∆ η ) . The termsarising when deriving the first and third-order differences are ( Λ Ti − − Λ Ti − )( κ i − / k i − / ) + ( Λ Ti − − Λ Ti )( k ν i − / κ i − / − k i + / κ i − / )+ ( Λ Ti − Λ Ti + )( − k ν i + / κ i + / + k i + / κ i + / ) + ( Λ Ti + − Λ Ti + )( − κ i + / k i + / ) , and may be rewritten as ( Λ Ti − − Λ Ti − ) κ i − / k i − / − ( Λ Ti − − Λ Ti ) κ i − / k i − / + ( Λ Ti − Λ Ti + ) κ i + / k i + / − ( Λ Ti + − Λ Ti + ) κ i + / k i + / − k (( Λ Ti − − Λ Ti ) ν i − / κ i − / − ( Λ Ti − Λ Ti + ) ν i + / κ i + / ) . The last two terms that stem from first-difference dissipation fluxes, are both O ( ∆ ξ ∆ η ) as the sensor in ξ ( I -mesh) direction is O ( ∆ ξ ) . The possibly problematic terms are the four first terms stemming from the third-differencedissipation fluxes. These are individually O ( ∆ ξ ∆ η ) terms and the question is whether their linear combination is ahigher order term. If the algebraic expression of k as function of the state variables is the same for all four values,then a third order difference is identified and linear combination of these four terms is O ( ∆ ξ ∆ η ) . If some of the ν l + / are equal to the left value ν l and some to the right value ν l + but none of the ¯ k is zero, the terms stemmingfrom the derivation of third-order differences may be rewritten ( Λ Ti − − Λ Ti − ) κ i − / k − ( Λ Ti − − Λ Ti ) κ i − / k + ( Λ Ti − Λ Ti + ) κ i + / k − ( Λ Ti + − Λ Ti + ) κ i + / k − ( Λ Ti − − Λ Ti − ) κ i − / k ν i − / + ( Λ Ti − − Λ Ti ) κ i − / k ν i − / − ( Λ Ti − Λ Ti + ) κ i + / k ν i + / + ( Λ Ti + − Λ Ti + ) κ i + / k ν i + / .
9n this case, the sum of the terms involving k is O ( ∆ ξ ∆ η ) and the terms in k may be individually identified as O ( ∆ ξ ∆ η ) terms. For unconsistency to be observed, it is hence necessary that: (a) due to the max in equation (20),the actual formulas of the involved k are not the same up to a shift of subscripts; (b) some of the k are zero. (Notethat these conditions are expected to appear at a marginal number of cells in particular as ν is very close to zeroexcept in the vicinity of shockwaves. Also note that these conditions are not event sufficient for dual unconsistency –if k i − / = k i + / = k i + / (cid:44) k i + / (cid:44) In a penultimate cell of subscript ( , j ) , the previous equations are slightly altered as fluxes F / , j and F / , j havespecific definitions involving a boundary state, W b / , j . Actually W b / , j is derived only from the adjacent conservativevariables W , j for each boundary interface; W b satisfies all Dirichlet like relations to be imposed at the boundary. Thenumerical fluxes at interfaces 3 / / ( , j ) (see Fig. 1), note first that the contributionsof the centered flux are standard for both mesh direction since F JST / , j does not depend on W , j . The terms of lines(26)-(27) are the same as in a current cell and do not need to be investigated again. The examination of those of lines(24)-(25) is done first considering the derivation of the sensor and the spectral radius. As the resulting contributionsto the discrete adjoint equation have been bounded individually in the previous subsection, it is only needed to lookat the specific contributions of faces ( / , j ) and ( / , j ) . No artificial dissipation is introduced at ( / , j ) , whilederiving the sensors in F JST / yields ( Λ T − Λ T ) ∂ ν / ∂ W κ / (cid:16) − ( W − W ) − [ k / > ] ( W − W + W ) (cid:17) which is O ( ∆ ξ ∆ η ) . The corresponding terms for F JST / and F JST / are at least O ( ∆ ξ ∆ η ) as justified in the previoussubsection. When deriving the spectral radius, the only specific term, depending on the boundary condition, is the onestemming from F JST / , ( Λ T − Λ T ) ∂ κ / ∂ W (cid:16) − k ν / ( W − W ) + k / ( W − W + W ) (cid:17) . The first part of this expression is O ( ∆ ξ ∆ η ) (a more precise assertion depends whether ν / = ν (that is O ( ∆ ξ ) ) or ν / = ν that is O ( ∆ ξ ) ). The second part of this expression is O ( ∆ ξ ∆ η ) . None of them hence prevents dualconsistency. Finally, the terms arising from the derivation of the differences terms in (24)-(25) read ( Λ T − Λ T )( − k ν / κ / − k / κ / ) + ( Λ T − Λ T )( k ν / κ / + k / κ / ) + ( Λ T − Λ T )( − κ / k / ) . They may be rewritten as − ( Λ T − Λ T ) k / κ / + ( Λ T − Λ T ) k / κ / − ( Λ T − Λ T ) κ / k / − ( Λ T − Λ T ) k ν / κ / + ( Λ T − Λ T ) k ν / κ / . The two-term difference of the last line is O ( ∆ ξ ∆ η ) . The sum of the terms of the first line is only O ( ∆ ξ ∆ η ) and introduces an inconsistency. Compared to the current-point equation, due to the specific F / , j formula, a ( Λ T − Λ T ) k / κ / is missing to have a third-difference appear.How the scheme can be corrected and the discrete adjoint can be slightly modified to obtain dual-consistency isdiscussed in § 4.4. 10 igure 1: Location of the state variable W and the corresponding boundary variable W b for the ( , j ) cell. Stencil of the JST scheme for this cell. Considering an output of the form (1) that depends on the wall boundary variable, W b , the discrete adjoint equationin a cell adjacent to a wall-boundary reads (see cell ( , j ) in Fig. 1) ∑ k Λ Tk ∂ R k ∂ W ( , j ) = − ∂ J ∂ W b ( / , j ) ∂ W b ( / , j ) ∂ W ( , j ) . (31)Due to the stencil of JST scheme in (19), (22) and (23), the k subscript varies in { ( , j − ) , ( , j − ) , ( , j ) , ( , j + ) , ( , j + ) , ( , j ) , ( , j ) } (see Fig. 1). The expanded expression of (31) reads ( Λ T ( , j − ) − Λ T ( , j − ) ) ∂ F , j − / ∂ W ( , j ) + ( Λ T ( , j − ) − Λ T ( , j ) ) ∂ F , j − / ∂ W ( , j ) (32) +( Λ T ( , j ) − Λ T ( , j + ) ) ∂ F , j + / ∂ W ( , j ) + ( Λ T ( , j + ) − Λ T ( , j + ) ) ∂ F , j + / ∂ W ( , j ) (33) +( Λ T ( , j ) − Λ T ( , j ) ) ∂ F / , j ∂ W ( , j ) + ( Λ T ( , j ) − Λ T ( , j ) ) ∂ F / , j ∂ W ( , j ) (34) − Λ T ( , j ) ∂ F / , j ∂ W ( j , ) = − ∂ J ∂ W b ( i , / ) ∂ W b ( i , / ) ∂ W ( , j ) . (35)It is clear that the terms of line (35) are O ( ∆ η ) . The J derivative fluxes appearing in (32) (33) are standard fluxesin the i -direction of the form (19) that do not depend on the boundary condition. The maximum order of magnitudeof the various terms of these lines is hence O ( ∆ ξ ∆ η ) (obtained when differentiating the centered flux). Whatever thespecific complexity of the ( ∂ F i , / / ∂ W ( , j ) ) terms due to the boundary state W b ( W ( , j ) ) , all terms appearing in line(34) are at least O ( ∆ ξ ∆ η ) . Keeping only the O ( ∆ η ) terms, previous equations yields − Λ T ( i , ) Sx ( / , j ) Sy ( / , j ) ∂ p b ( / , j ) ∂ W b ( / , j ) ∂ W b ( / , j ) ∂ W ( , j ) (cid:39) − ∂ (( − S ( / , j ) · d ) p b ( / , j ) ) ∂ W b ( / , j ) ∂ W b ( / , j ) ∂ W ( , j ) . That is simplified in Λ ( i , ) Sx ( / , j ) + Λ ( i , ) Sy ( / , j ) (cid:39) − S ( / , j ) · d . (36)Equation (36) is obviously the cell-centered FV counterpart of the continuous adjoint wall boundary condition (8)in the Cartesian frame. Note that it involves metric terms at the wall and the adjoint solution at adjacent cells and alsothat it has been derived by equalizing lower-order (i.e., first-order) terms in space. How accurately equation (36) issatisfied is discussed in sections 6.1 to 6.3. 11 .4. Modified discrete adjoint for adjoint consistency The discrete adjoint of the cell-centered JST scheme (19), (22) and (23) is thus adjoint-consistent except at penul-timate cells w.r.t. a boundary. This issue is considered here. A satisfactory modification of F JST / for dual consistency,would be F / = ( F ( W ) + F ( W )) · S / − k ν / κ / ( W − W ) + k / κ / ( W − W + W ) Unfortunately, as there is no artificial dissipation in F / so that the k / κ / ( W − W + W ) term of this modified F / would introduce an anti-dissipative term in the equivalent equation of the cell adjacent to the boundary that mayresult in a significant loss of robustness. This variant has hence not been considered.Adjoint-consistency can otherwise be recovered at the price of a minor and local approximation in the Jacobianof the scheme. Following the calculations of § 4.2, the needed change in the differentiation of the F / , j flux is thereplacement of − ( Λ T − Λ T ) k / κ / by − ( Λ T − Λ T ) k / κ / in ( ∂ F / , j / ∂ W ) when differentiating the third order difference. Although all adjoint values are coupled, this induces no significantchange far from the boundaries. Close to the wall, this change makes the contours more regular and lowers themaximum values of the adjoint field as observed in Fig. 2 for a transonic flow about the NACA0012. Besides itis noted that adjoint contours are stretched close de penultimate cells w.r.t. the boundary for the inconsistent exactdiscrete adjoint (Fig. 2 left). Moreover, when plotting the residual of the continuous adjoint equation discretized onthe mesh with the computed flow and discrete discrete adjoint (see Fig. 3 for the same transonic flow), lower valuesare observed for the consistent adjoint variant.Due to these observations, the slightly modified and dual-consistent discrete adjoint is used when discussing meshconvergence and asymptotic behavior at the stagnation streamline and at the wall. In all other sections, the exactdiscrete adjoint solution of the scheme is retained. Figure 2: ( M ∞ = . α = o , 2049 × Λ CL p : inconsistent discrete adjoint (left) and consistent discrete adjoint (right).
5. Adjoint Rankine-Hugoniot relations
We consider an isolated discontinuity Σ where the direct solution satisfies the RH jump relations [[ n x F x ( W ) + n y F y ( W )]] = , (37)12 igure 3: ( M ∞ = . α = o , 2049 × where [[ · ]] denotes the jump across Σ and n = ( n x , n y ) T is the normal to Σ . Introducing t = ( − n y , n x ) T , the adjointequation (4) may be rewritten as − A Tn ∂ ψ∂ n − A Tt ∂ ψ∂ t = Ω , (38)with A n = n x A + n y B = n x n y ( γ − ) e c n x − uv n v n − ( γ − ) un x un y − ( γ − ) vn x ( γ − ) n x ( γ − ) e c n y − vv n vn x − ( γ − ) un y v n − ( γ − ) vn y ( γ − ) n y (cid:0) ( γ − ) e c − H ) v n Hn x − ( γ − ) uv n Hn y − ( γ − ) vv n γ v n , A t = − n y A + n x B = − n y n x − ( γ − ) e c n y − uv t v t + ( γ − ) un y un x + ( γ − ) vn y − ( γ − ) n y ( γ − ) e c n x − vv t − vn y + ( γ − ) un x v t − ( γ − ) vn x ( γ − ) n x (cid:0) ( γ − ) e c − H ) v t − Hn y − ( γ − ) uv t Hn x − ( γ − ) vv t γ v t , where v n = n · U , v t = t · U , ∂ n ψ = n x ∂ x ψ + n y ∂ y ψ , and ∂ t ψ = − n y ∂ x ψ + n x ∂ y ψ . Assuming that A n is non-singular,the adjoint is continuous across Σ [36, 1, 38, 37]: [[ ψ ]] =
0, so [[ ∂ t ψ ]] =
0. Likewise we have [[ ρ v n ]] = [[ v t ]] =
0, and [[ H ]] = + u x (38b) + u y (38c) gives v n (cid:0) ∂ n ψ − H ∂ n ψ (cid:1) + v t (cid:0) ∂ t ψ − H ∂ t ψ (cid:1) = , and applying the jump operator we obtain [[ v n ∂ n ψ ]] − H [[ v n ∂ n ψ ]] = . (39)(For functions of the static pressure and the geometry, according to equation (16), we expect the left and right terms ofthe difference to be zero which is actually observed later on in the numerical experiments with lift and drag adjoint.)Likewise, (38d) gives ( γ − )[[ n x ∂ n ψ + n y ∂ n ψ ]] + γ [[ v n ∂ n ψ ]] = . (40)13ow, the operations [[ n x (38c) − n y (38b) ]] and [[ n x (38b) + n y (38c) ]] provide [[ v n (cid:0) − n y ∂ n ψ + n x ∂ n ψ + v t ∂ n ψ (cid:1) ]] + [[ v n ]] (cid:0) n x ∂ t ψ + n y ∂ t ψ (cid:1) = , (41a) [[ ∂ n ψ ]] + [[( u + v n n x ) ∂ n ψ + ( v + v n n y ) ∂ n ψ + ( H + v n ) ∂ n ψ ]] + [[ v n ]] v t ∂ t ψ = . (41b)Finally, the adjoint equation to the RH relations (37) reads [[ − n y F x ( W ) + n x F y ( W )]] · ∂ t ψ = v t [[ ρ ]]( ∂ t ψ + H ∂ t ψ ) + ([[ p ]] + v t [[ ρ ]])( − n y ∂ t ψ + n x ∂ t ψ ) = . (42)
6. Numerical experiments
We consider the direct and adjoint solutions for inviscid flows around the NACA0012 airfoil in different regimes.Steady-state solutions are computed on a series of five structured meshes (ranging from 129 ×
129 to 4097 × elsA code [55]. The farfield boundary is set at about 150 chords. Thecomputations have been converged with a second order accurate finite volume method using the JST scheme (19),(22) and (23) with k = . k = . k = k = . elsA [57]for the two functions and the complete set of meshes. We first consider a supersonic flow about the NACA0012 airfoil with a free-stream Mach number of M ∞ = . α = o .The solution is displayed in Fig. 4. The flow is supersonic and constant up to a detached shock wave. Downstreamthe shock wave, the flow is subsonic in a small bubble close to the airfoil leading edge and supersonic elsewhere. Itaccelerates along the airfoil up to a fishtail shock wave based on the trailing edge. Downstream this second shockwave, the flow is still supersonic with a Mach number close to the upwind farfield Mach number. Figure 4: ( M ∞ = . α = o , 2049 × δ R source term ( ε =6 . − )at point ( x = . , y = . ) . .1.1. Lift and drag adjoint solutions Figure A.26 presents the contours of all four components of the adjoint fields.The gradient of the discrete adjoint is then estimated at cell centers using Green formula and finally the residualin cell ( i , j ) of the continuous adjoint equation, res i j = − A Ti j (cid:18) ∂ Λ ∂ x (cid:19) i j − B Ti j (cid:18) ∂ Λ ∂ y (cid:19) i j , (43)is evaluated for both functions. In the above equation, A i j and B i j are the Jacobians evaluated with cell-centervalues. For a first verification, both terms of the residual and their sum are plotted separately for the lift, on the2049 × | res | below a threshold is decreasing as the mesh is refined. Nevertheless significant (andincreasing) | res | values are still observed on the finest meshes at the trailing edge and along the Mach lines passingthrough the trailing edge (where a perfect adjoint field would exhibit a discontinuity due to the end of the functionsupport) – see Fig. 6.Concerning farfield boundary conditions, the normal Mach number appears to be supersonic at the intersection ofthe characteristic geometrical strips (see Fig. 7) and the far field boundary. No continuous adjoint boundary conditionsis hence to be applied there and no check of discrete versus continuous adjoint is required. Concerning the discreteequation (36) that is the counterpart of the continuous boundary condition at the wall , equation (8), its two terms areplotted in the right part of Fig. 6. They appear to be superimposed so that (36) is actually well satisfied.Besides, as a mean to check the mesh convergence of adjoint fields, we have calculated the finest mesh adjoint-lift and adjoint-drag on another 2049 × Figure 5: ( M ∞ = . α = o , 2049 × − A T ( ∂ Λ CL p / ∂ x ) discretized with discrete adjoint field (first component of the vector).Right: − B T ( ∂ Λ CL p / ∂ y ) discretized with discrete adjoint field (idem). The flow is supersonic and constant upstream of the detached shock wave so that the Jacobians A and B in thecontinuous adjoint equation (4) are constant matrices in this area that enables a specific analysis of the continuousadjoint solution in this region. This analysis, partly presented in [42] is here further developed and assessed. In asupersonic regime, the direct and adjoint equations both exhibit specific directions of propagation corresponding tosimple waves solutions to (4): ψ ( x , y ) = φ ( x sin ( γ ) − y cos ( γ )) λ , where γ is the angle made by the direction of propagation with the x-axis, λ is a vector representing the convectedinformation and φ is a scalar function. Injecting this expression into (4) yields15 igure 6: ( M ∞ = . α = o ) Left: Residual of the continuous equation (first component) evaluated with the lift discrete adjoint fields (bottom:129 ×
129 mesh, top: 2049 × × φ (cid:48) ( x sin ( γ ) − y cos ( γ )) × ( sin ( γ ) A T − cos ( γ ) B T ) λ = . This equation admits a non-trivial solution λ if and only ifdet (cid:0) sin ( γ ) A T − cos ( γ ) B T (cid:1) = . This condition is the same as the for supersonic Euler equations since the transposition plays no role in thecalculation of the determinant. The eigenvalues of sin ( γ ) A − cos ( γ ) B are: ( u x sin ( γ ) − u y cos ( γ ) , u x sin ( γ ) − u y cos ( γ )+ c , u x sin ( γ ) − u y cos ( γ ) − c ) . Since one of the eigenvalues vanishes, γ is solution of one of these equations: u x sin ( γ ) − u y cos ( γ ) = , u x sin ( γ ) − u y cos ( γ ) + c = , u x sin ( γ ) − u y cos ( γ ) − c = . Let α denote the angle of incidence, then the velocity components then read u x = (cid:107) U (cid:107) cos ( α ) and u y = (cid:107) U (cid:107) sin ( α ) .Substituting these expressions in the above relations yields γ = α or α + πγ = α − µ or γ = α − µ + πγ = α + µ or γ = α + µ + π , where µ = arcsin ( / M ∞ ) is the Mach angle. Hence, the information propagates along the three privileged directions γ = α , γ = α − µ et γ = α + µ . Evidently, the domains of influence/dependance are inverse one another for state andadjoint variables.The λ eigenvectors of interest are those of sin ( γ ) A T − cos ( γ ) B T . These are also the left eigenvectors that appearin the more classical diagonalization of the inviscid flux Jacobian. Following [49], the (left) eigenvector associated to γ = α − µ (null eigenvalue u x sin ( γ ) − u y cos ( γ ) + c ) and γ = α + µ (null eigenvalue u x sin ( γ ) − u y cos ( γ ) − c ) are λ α − µ = c ρ (cid:16) − u x n x − u y n y c + ( γ − ) M (cid:17) ρ (cid:16) n x c − ( γ − ) u x ρ (cid:17) ρ (cid:16) n y c − ( γ − ) u y ρ (cid:17) γ − ρ c , λ α + µ = c ρ (cid:16) u x n x + u y n y c + ( γ − ) M (cid:17) ρ (cid:0) n x + ( γ − ) u x c (cid:1) ρ (cid:0) n y + ( γ − ) u y c (cid:1) γ − ρ c , ( n x , n y ) T , which may be replaced by ( sin ( γ ) , − cos ( γ )) to recover the solution as a function ofthe simple wave angles. Checking that the simple waves solution satisfy relation (16) uses the nullity of the eigenvalue( u x n x + u y n y + c = u x n x + u y n y − c = γ = α is two but the classically exhibited left eigenvectors [49] donot satisfy the Giles and Pierce relation for pressure-based outputs (16). The vectors of this 2D space that satisfy thisrelation is the 1D space generated by λ α = − − ( γ − ) q c ( γ − ) u x c + n y n y u x − n x u y ( γ − ) u y c − n x n y u x − n x u y − γ − ρ c . Figure 7 highlights the directions of propagation in the first component discrete adjoint field.The ratio of theadjoint components inside all three geometrical strips is found equal to the ratio of the corresponding componentsof the associated λ vector. The final demonstration of the consistency between the discrete adjoint fields and thetheoretical form of the continuous adjoint solution upwind the shock wave, is obtained by performing a line extractionof the three terms in ψ ( x , y ) = φ α ( x sin ( α ) − y cos ( α )) λ α + φ α − µ ( x sin ( α − µ ) − y cos ( α − µ )) λ α − µ + φ α + µ ( x sin ( α + µ ) − y cos ( α + µ )) λ α + µ , (44)and propagating these values in the complete zone upwind the detached shock wave. More precisely the values ofthese terms are extracted from a fine discrete adjoint field at a distance 2 . c from the leading edge (where the threebands are separated). They are then interpolated everywhere upwind the shock wave according the local values of x sin ( α ) − y cos ( α ) , x sin ( α − µ ) − y cos ( α − µ ) and x sin ( α + µ ) − y cos ( α + µ ) . For the finer meshes, the resultingfield appears to be visually identical to the actual discrete adjoint field, just upstream of the shock wave, where thebands are superimposed (results not shown here). Figure 7: ( M ∞ = . α = o , 4097 × CL p adjoint. .1.3. Residual perturbation analysis If the residual perturbation is located in a supersonic zone, its influence on the flow close to the source appearsto be restricted to the downstream part of the streamline and of the two characteristic curves passing by the sourcelocation. Moreover, assuming a constant flow W , by the same arguments as in previous section, its perturbation δ W can be expressed as δ W ( x , y ) = ψ α ( x sin ( α ) − y cos ( α )) r α + ψ α − µ ( x sin ( α − µ ) − y cos ( α − µ )) r α − µ + ψ α + µ ( x sin ( α + µ ) − y cos ( α + µ )) r α + µ , (45)where ψ γ stands for a small perturbation scalar function and each r γ vector is a right eigenvector of the Jacobiansin ( γ ) A − cos ( γ ) B . The consistency of the actual numerical δ W with this analytical model can be simply assessed,upwind the detached shock, for δ R and δ R , since they activate the Mach-lines part of the solution, ψ α − µ ( x sin ( α − µ ) − y cos ( α − µ )) r α − µ + ψ α + µ ( x sin ( α + µ ) − y cos ( α + µ )) r α + µ , with one-dimensional vector spaces span { r α − µ } and span { r α − µ } : it was successfully verified that the ratio of thecomponents of the numerical δ W and δ W is the one of the corresponding right eigenvector (results not shownhere).On the contrary, if the source term is located in the subsonic bubble or if one of the three curves enters the subsonicbubble, then the support of flow perturbation includes this complete subsonic area consistently with the elliptic natureof the equations [58, chap. 11].More details are given for the case when the source term is located in the supersonic zone. Let consider theperturbation decomposition (45). If the entropy or the total enthalpy (that are convected in steady state Euler flows)are perturbed ( δ R or δ R ), close to the source, the perturbation is located along the trajectory downstream the sourceand possibly along shock waves for δ R . If these quantities are not altered by the perturbation ( δ R or δ R ), theresponse of the steady state flow to the perturbation, close to the source, is located along the two characteristic linesdownwind the source and the usual relative-variations along trajectories of thermodynamic and kinetic variables arevalid between nominal and locally perturbed flow. More precisely : • δ R induces no change in the stagnation quantities and entropy. Mach number is decreased along the two char-acteristic lines starting from the perturbation location. The static pressure, temperature and density are increasedand velocity components are also perturbed along these two lines with a decrease of velocity magnitude; • δ R induces no change in the stagnation quantities and entropy. Along the upper characteristic curve startingfrom the perturbation point, the static pressure, density, temperature are increased whereas Mach number, ve-locity magnitude, x-component of velocity are decreased. Opposite variations are observed along the lowercharacteristic curve. The y-component of the velocity is increased along both characteristic curves; • δ R induces no variation of stagnation pressure, static pressure and Mach number. An increase of total enthalpy(stagnation temperature), temperature, both components of velocity and entropy is observed along the trajectorystarting from the perturbation location. The density is decreased; • δ R induces no local variation of the static pressure except along shock waves and no variation of the totalenthalpy. Mach number, density, stagnation pressure, stagnation density, and both velocity components areincreased along the trajectory starting from the perturbation location. Entropy and temperature are decreasedalong the trajectory. As mentioned end of § 3.3 this notion of perturbation downstream the source and along astreamtube of the initial is only a main order description as the streamtubes structure of a non-constant flow isperturbed by a δ R source. all variations are given for a positive ε in equation (13) and positive components of velocity for the nominal flow δ W perturbation on the pressure fieldplotted together with the characteristic curves passing by the source location. The adjoint variables are continuous across shocks, whereas their derivatives may be discontinuous [36, 1, 38].In Fig. 8, a clear discontinuity of the adjoint-drag or the adjoint-lift gradient is observed for the adjoint componentassociated to z -coordinate momentum equation. Figure 8: ( M ∞ = . α = o , 4097 × z -coordinate momentum equation. The positionof the shock is indicated by a pink line. Figure 9 validates the adjoint RH jump relations (39) to (41b) derived in § 5. We here display both lift and dragadjoint solutions. The shock location is about s = .
11 and we observe that the adjoint RH relations are well satisfied,while the normal gradients of the adjoint variables display large discontinuities at the shock position.
We now consider a transonic flow around the NACA0012 airfoil with conditions M ∞ = .
85 and α = o . A strongshock wave develops on the suction side and a weaker shock on the pressure side (see Fig. 10). Figure A.27 presents the contours of all four components of the drag and lift adjoint fields. (As indicated before,the choice of the scale for the contours of first and last adjoint components in Figs. A.26, A.27 and A.28 allows tovisually check that equation (16) is closely satisfied by the discrete adjoint fields.) Both fields exhibit strong valuesand gradients close to the wall, the stagnation streamline, the two characteristics impinging the upper side and lowershock foot and two other characteristic curves that form a circumflex with the previous two. These features have beenobserved previously [8, 42, 39, 38].The physical point of view of Giles and Pierce [1] recalled in § 3.3 is again considered and the δ CL dp termscorresponding to the four source terms in (13) are plotted in Figs. 10 (right) and A.20. We observe that δ R and δ R have a strong impact on CL p when located along the characteristic curves that form the two circumflex. δ R has astrong impact when it is located close to the wall or to the stagnation streamline. The asset of this approach is obviouswhen studying the singular behavior of the adjoint fields: the δ J d responses to the physical source terms exhibit onlypart of the singular zones of the classical adjoint fields and the high values/gradients of one δ J d are transfered to allor to several of the usual adjoint components according to equation (15).The residual of the continuous adjoint equation is calculated as detailed in § 6.1.1. As the mesh is refined, thezones with significant residual have a decreasing area but remain large. Increasing values are also observed close to19 5 · − . .
15 0 . . .
52 shock s MH ρ v n v t · − . .
15 0 . − s ∂ n ψ ∂ n ψ ∂ n ψ ∂ n ψ · − . .
15 0 . − s (39)(40)(41a)(41b) · − . .
15 0 . − − s (39)(40)(41a)(41b)Figure 9: ( M ∞ = . α = o , 4097 × ( − . , . ) and ( . , . ) as a function of the local coordinate s . An equation number refers tothe arguments between brackets [[ · ]] in the corresponding equation, e.g., v n ( ∂ n ψ − H ∂ n ψ ) has been plotted for (39).Figure 10: ( M ∞ = . α = o , 2049 × δ CL p (equations (13) and (14) with ε =1),location of the points for the residual perturbation analysis plus sonic line (in pink). λ terms are most probably too large for a simple main-order in space analysis. Figure 11: ( M ∞ = . α = o ) Left: residual of the continuous equation (first component) evaluated with the lift discrete adjoint fields (bottom:129 ×
129 mesh, top: 2049 × × Figure 12 highlights a clear discontinuity across the shock of the gradient of the drag or the lift adjoint componentassociated to z -coordinate momentum. This is further analyzed in Fig. 13 where the lift and drag adjoint derivativesare plotted along a line crossing the shock. Again it is observed that the adjoint RH relations (39) to (41b) are wellsatisfied, while the normal adjoint derivatives may be discontinuous across the shock. δ R and δ R perturbations along specific characteristic lines The first order variation of CL p and CD p in response to the physical δ R perturbations in (13) and (11) are calculatedand plotted for CL p in Figs. 10 (right) and A.20. More precisely, we set ε = δ J iscalculated through (11) and affected in the plot to the source location.The δ R and δ R sources, lead to δ CLp and δ CLp δ CDp and δ CDp ) exhibiting strong (positiveand negative) values along characteristic lines of the supersonic areas upwind both shocks, in particular along thecharacteristic lines impacting the shock foots. For a better understanding, the flow perturbation and the lift perturbation δ CLp due to the δ R source terms are calculated for sources located in the upper side supersonic bubble at points 4,5, 6 and 7 (see Fig. 10). The changes in the wall pressure field (of the re-converged flow) due to δ R source termsare illustrated in Fig. 14. From the discussion in § 6.1.3, it is known that a positive pressure increment is propagatedalong the two Mach lines starting from the source location. Its influence at the wall is observed along with the one ofother more global perturbations in the flow propagated through the subsonic zones: for a source term located at point4, there is a visible local increase (and oscillation) of the pressure at about x = .
11 where the lower characteristicpassing by point 4 reaches the wall. Nevertheless, the change in
CLp is mainly due to global variations of the staticpressure all along the airfoil. For point 6, the change in
CLp is due to the displacement of the lower-side shock andthe variation in pressure where the lower characteristic passing at point 6 impacts the wall (small bump at x (cid:39) CLp is due to the displacement of both, the upper side andthe lower side shock. Point 5 has been placed close to point 6 and 7 in order to better understand why a point closeto the two well-marked characteristic curves may have a weak influence on
CLp and
CDp . Actually, for this point,the influence of the upper side shock displacement and the lower characteristic impact (small bump at x (cid:39) δ CLp value.21 igure 12: ( M ∞ = . α = o , 4097 × z -coordinate momentum equation (sonic line inpink). · − . .
15 0 . s MH ρ v n v t · − . .
15 0 . − s ∂ n ψ ∂ n ψ ∂ n ψ ∂ n ψ · − . .
15 0 . s (39)(40)(41a)(41b) · − . .
15 0 . s (39)(40)(41a)(41b)Figure 13: ( M ∞ = . α = o , 4097 × ( . , . ) and ( . , . ) as a function of the local coordinate s . An equation number refers to thearguments between brackets [[ · ]] in the corresponding equation, e.g., v n ( ∂ n ψ − H ∂ n ψ ) has been plotted for (39). igure 14: ( M ∞ = . α = o , 2049 × δ R (equation (13) with ε = . e − ) source termslocated at points 5, 6, 7 and 8 (variation w.r.t. nominal flow 50 times amplified) δ R perturbation in the vicinity of the wall and the stagnation streamline As mentioned beforehand observing Fig. A.20, δ R is the physical source term responsible for high values oflift/drag-adjoint and gradient of adjoint close to the wall and the stagnation streamline. The mechanism by which δ R modifies the flow field and the forces values is detailed for the vicinity of the stagnation streamline and the wall upperside.In the first series of considered perturbed flows, the δ R source terms are located upper side, in the supersonicarea, along a mesh line orthogonal to the wall, close to x = .
5. In a supersonic zone, the perturbation of the flowcaused by a δ R source term along the corresponding streamline has been described in section § 6.1.3. Actually, atlarge scale, the pressure field is not altered except by the displacement of the two shocks and the upper side shockposition of the nominal flow is not compatible with the increased upwind Mach number (ie for this almost normalshock, the classical equation for the ratio of downwind over upwind static pressure as a function of the upwind Machnumber can not be satisfied where the perturbed streamline hits the shock). In fact, the upper side shock is movedglobally backwards and the shifted shock goes with a decreased upwind static pressure and an increased downwindstatic pressure (see local slopes of the orange curve in Fig. 14). It is observed that the smaller the distance of thesource and perturbed streamline to the wall, the larger the displacement of the upper side shock foot which explainsthe growth of δ CLp in the vicinity of the wall. Observing the static pressure changes at smaller scales, it is noted,as expected, that the perturbation of the static pressure in the subsonic area is not restricted to the continuation of thestreamline but extends to the whole subsonic area close to the profile downstream the shocks. As concerning the staticpressure close to the lower side, at these finer scales, an increase is observed upwind the shock and below the trailingedge whereas a decrease is observed downwind the shock; this goes with a backward displacement of the lower shock.This is illustrated in Fig. 15 (left). Both changes of the static pressure at the wall correspond to an increase of the liftthat is consistent with the local values of δ CLp .If the δ R source is located close to the wall but downwind the upper side shock, similar shock displacementsand lower side fine scale perturbations of the pressure are observed. On the contrary, the flow perturbations createddownwind the upper side shock are very different for sources located on the same streamline either sides of the shock.Concerning the points located close to the stagnation streamline, the mechanism that produces the lift perturbation δ CLp due to the δ R source, seems to be correlated with particles convection since the iso- δ CLp roughly follow thestreamlines (see Fig. A.20 right). Independently of the equations specific to supersonic flows or supersonic bubbles in23 igure 15: ( M ∞ = . α = o , 2049 × δ R source terms located close to the wall on shocks position (source atpoints ( , j ) j = , , , , ,
35) – right: Influence of δ R source terms located close to the stagnation streamline on shocks position (source atpoints ( + l , ) l = − , − , − , − , , , , , transonic flows that have been used before, let us first note that the base transonic flow satisfies everywhere d p / dt = δ R is obtained by the partial derivation of theinviscid flux in the flow direction with respect to the stagnation pressure, at constant total enthalpy and static pressureand, from a fluid dynamics point of view, the expected main order perturbation (see end of §3.3) of the stagnationpressure field is the convection of the local increment created by the source downwind the source location. For anumerical point of view, a dissipated propagation of the δ p from the source location is observed (see Fig. A.24that also illustrates the propagation of the decrease of entropy created by the source). After this increment has beenpropagated up to the upper side shock, the mechanism through which both shocks are moved is the same as the onedescribed before. The perturbation of the pressure distribution is illustrated in Fig. 15 (right) for a series of sourceslocated close to the stagnation streamline, about 0 . c upwind the stagnation point.In [1], the authors consider a 2D inviscid flow about an airfoil and an output J defined as the integral along thewall of the pressure times a local factor. With the additional approximation of potential flow, they then derive theasymptotic behavior of δ J m , the variation of J under δ R , when the source term is located in the vicinity of thestagnation streamline. Subject to this assumption, they prove that δ J (cid:39) ± d − / where d is the distance of the source m to the stagnation streamline and the sign is changed when crossing the streamline. This relation is not alwaysnumerically well satisfied by numerical lift / drag adjoint fields of Euler flows [39, 38]. For the transonic test-caseof interest, it as been demonstrated above that δ CLp m and δ CDp m result in similar perturbations of the flow forpoints m located on the same streamline, close to the stagnation streamline and close to the wall upwind the shock.As the considered meshes are structured and essentially orthogonal at each vertex (in particular at the wall) [54], iteasy to search for the γ exponent such that δ J (cid:39) d γ in the vicinity of the wall. For our numerical solutions, thisexponent appears to depend on the x location (rather than on CD vs CL , lower side vs upper side or mesh density).It is decreasing with increasing x , from 0. at the leading edge to about -0.1 at x = . V ∧ ( ∇ × V ) = − T ∇ s and results in the following perturbation equation δ ( V ) ∧ ( ∇ × V ) + V ∧ δ ( ∇ × V ) = − δ ( T ∇ s ) , (46)that appears to be well satisfied by numerical solutions (see e.g. Fig. A.25). As recent articles discussed the singularityof lift and drag adjoint at the trailing edge, as δ R is the physical source term responsible for high values / high gradientof adjoint at the wall, we finally note that the perturbed Crocco equation is singular at the trailing edge (and moregenerally along the trailing edge slip line): it is the only location were the δ ( V ) ∧ ( ∇ × V ) term is not negligible withrespect to V ∧ δ ( ∇ × V ) which may enter the discussion about the singularity of lift/drag adjoint at the trailing edge.24 .3. Subsonic regime Finally, we consider a smooth subsonic flow with freestream conditions M ∞ = . α = o . The Mach numbercontours are presented in Fig. 16 where we observe the stagnation point located just under the leading edge as well asa strong acceleration of the flow on the upper side due to the profile incidence. Figure 16: ( M ∞ = . α = o , 2049 × Figure A.28 presents the contours of all four components of the adjoint fields. The lift adjoint components exhibitstrong values and gradients close to the wall and the stagnation streamline. The amplitudes of these values increaseas the mesh is refined. Conversely, all components of the drag adjoint appear to be bounded and mesh converging.These different behaviors may be attributed to the facts that the adjoint field quantifies the sensitivity of the quantityof interest w.r.t. the residuals, and that the limit inviscid drag is zero for a subsonic flow whereas the limit lift is not.The residuals of the continuous adjoint equations have been calculated as detailed in § 6.1.1. As the mesh isrefined, the zones with significant residuals for the lift-adjoint have a decreasing support but exhibit increasing valuesclose to the wall and to the stagnation streamline (see left part of Fig. 17). The farfield adjoint boundary condition issatisfied by the discrete adjoint fields as the adjoint field is almost zero at the farfield boundary. Right part of Fig. 17shows that the adjoint wall boundary condition (36) is well satisfied for both lift and drag.The physical point of view of Giles and Pierce [1] is then considered to gain understanding for the behavior ofthe lift-adjoint . The δ CLp , δ CLp , δ CLp , and δ CLp responses have been calculated and plotted in Fig. A.21.We observe that δ R is the only source responsible for the singular behavior of the adjoint-lift close to the wall andstagnation streamline. As for transonic or supersonic flows, a δ R source results in positive increments of the stagnation pressure andthe stagnation density and a negative increment of entropy, and all three increments are convected downstream thesource location following the classical laws for Euler flows with uniform total enthalpy. However, the response of the let us briefly note that this approach is somehow making the analysis more complex for this flow condition when applied to the drag. Actually,the numerical solutions exhibit dCD values such that dCD → − ∞ in the vicinity of the stagnation point (see right subplot of Fig. A.22) but notethat all terms of the last line of the matrix in equation (15) vanish. This finally results in bounded and mesh converging adjoint components for drag igure 17: ( M ∞ = . α = o ) Left: residual of the continuous equation (first component) evaluated with the lift discrete adjoint fields (down129 ×
129 mesh, up 2049 × × subsonic flow to the local perturbation of these three quantities, is very different from the one described in § 6.2.4for the selected transonic flow: density, velocity, static pressure, Mach number and all derived variables exhibit aglobal perturbation (see Fig. A.24). When the source is located close to stagnation streamline, the perturbation is thestrongest close to the leading edge. If it is located just under (resp. above) the stagnation streamline, the main effecton the static pressure at the wall is a decrease all along the pressure (resp. suction) side consistently with the localsign of δ CLp (Fig. A.21 right and Fig. 18).As in the transonic case, we examine whether δ CLp (cid:39) d γ in the vicinity of the wall. Once again, a negative γ exponent decreasing from 0 to about -0.20 with increasing x is derived from this assumption – see Fig. A.23 right. Figure 18: ( M ∞ = . α = o , 2049 × δ R source term located at points A, B, C, D, on the staticpressure at the wall. Distance of points A,B, C, D to the leading edge (cid:39)
7. Conclusion
In this work, we address open issues in the field of discrete and continuous adjoint equations for inviscid flows.The results are illustrated through in-depth numerical experiments on lift- and drag-adjoint fields for inviscid flows26bout the NACA0012 airfoil in supersonic, transonic and subsonic regimes.We first investigate the dual consistency of the JST scheme in cell-centered FV formulation for the discretization ofthe compressible Euler equations in 2D. We prove dual consistency of the scheme at interior cells, while inconsistencyoccurs at the penultimate cell to a physical boundary due to the absence of artifical dissipation at the boundary. Wehence propose a slight modification of the exact scheme derivation w.r.t. the penultimate-cell flow variables to recoverconsistency. The implementation of this modification is straightforward and was shown to improve the convergenceof the adjoint computation at wall boundaries under grid refinement. The modification has been carried out in ourcode and successfully illustrated. The link between continuous inviscid adjoint boundary condition at the wall andcell-centered discrete adjoint equation in the cell adjacent to the wall has also been established.Second, a new heuristic method has been proposed to discuss the adjoint consistency of discrete adjoint fields.It consists in discretizing the continuous equation with the discrete flow and discrete-adjoint field and discuss theconvergence of the residual towards zero as the mesh is refined. Coherently with the demonstrated property of adjointconsistency, this residual has been found to be vanishing except in areas where: (a) a discontinuous adjoint field isexpected (typically, vicinity of the trailing edge for supersonic flows); (b) the discrete adjoint field exhibits increas-ing values values and gradients and where theoretical divergence is suspected (stagnation streamline, wall, specificcharacteristic lines in supersonic areas).Third, we derive the general adjoint Rankine-Hugoniot relations across a shock and give relations linking thenormal derivatives across the shock. These relations have been numerically validated for a supersonic and a transonicflow. Besides, as the studied shockwave of the transonic flow is almost normal, the presented results answer questionsrecently raised about the adjoint gradient at shockwaves [38].We then use these tools for the analysis of the adjoint fields and in particular the analysis of the zones whereincreasing values and gradients of lift/drag adjoint are observed about profiles, at some flow conditions, when the meshis refined [1, 39, 38, 59]. The novelty in this domain is both, to take a detailed benefit of the source term approachof Giles and Pierce [1] identifying the δ R source term(s) associated to the known zones of numerical divergence, andto consider some corresponding shifted flows W + δ W (resulting of the δ R source) to understand the fluid dynamicsmechanisms that provoke strong changes in lift and drag. The focus has been put on the δ R source that appears tobe associated with the high values and gradients of the discrete adjoints, that are observed at many transonic flowconditions, close to the wall and to the stagnation streamline. This δ R source actually results in the convectionof an increment of stagnation pressure, stagnation density and entropy (provoking shock displacement) at globallyconstant total enthalpy and locally constant static pressure. Many properties of transonic flows lift or drag adjointfields are clarified by this approach: it first unifies the behavior expected at the wall and stagnation streamline. Italso explains the difference of numerical behavior in lift or drag adjoint at high transonic Mach numbers (shock-wave feet at trailing edge, no possible shock-foot displacement, no adjoint numerical divergence) and lower transonicMach numbers (more upwind shock-wave feet locations, shock feet displacement under the influence of δ R , adjointnumerical divergence). Besides, this approach provides a mechanical counterpart to the mathematical study of the liftand drag adjoint asymptotic behavior in the vicinity of the wall and stagnation streamline at the flow conditions forwhich numerical divergence is observed.Future work will concern the extension of this analysis to laminar flows considering the usual FV cell-centereddiscretizations of the viscous flux. Acknowledgments
The authors express their warm gratitude to J.C. Vassberg and A. Jameson for allowing the co-workers of D.Destarac to use their hierarchy of O-grids around the NACA0012 airfoil.This research did not receive any specific grant from funding agencies in the public, commercial or not-for-profitsector 27 iblographyReferences [1] Giles, M. and Pierce, N. Adjoint equations in CFD: Duality, boundary conditions and solution behaviour. In
AIAA Paper Series, Paper97-1850 . (1997).[2] Frank, P. and Shubin, G. A comparison of optimisation based approaches for a model computational aerodynamics design problem.
Journalof Computational Physics , 74–89 (1992).[3] Pironneau, O. On optimum design in fluid mechanics. J. Fluid Mech. (1), 97–110 (1974).[4] Jameson, A. Aerodynamic design via control theory. Journal of Scientific Computing (3), 233–260 (1988).[5] Jameson, A., Martinelli, L., and Pierce, N. Optimum aerodynamic design using the Navier-Stokes equations. Theoretical and ComputationalFluid Dynamics (1), 213–237 (1998).[6] Castro, C., Lozano, C., Palacios, F., and ZuaZua, E. Adjoint approach to viscous aerodynamic design on unstructured grids. AIAA Journal (9), 2125–2139 (2007).[7] Becker, R. and Rannacher, R. An optimal control approach to a posteriori error estimation in finite element methods. Acta Numerica ,1–102 (2001).[8] Venditti, D. and Darmofal, D. Grid adaptation for functional outputs: Application to two-dimensional inviscid flows. Journal of Computa-tional Physics , 40–69 (2002).[9] Dwight, R. Heuristic a posteriori estimation of error due to dissipation in finite volume schemes and application to mesh adaptation.
Journalof Computational Physics , 2845–2863 (2008).[10] Loseille, A., Dervieux, A., and Alauzet, F. Fully anisotropic mesh adaptation for 3D steady Euler equations.
Journal of ComputationalPhysics , 2866–2897 (2010).[11] Fidkowski, K. and Roe, P. An entropy approach to mesh refinement.
SIAM Journal of Scientific Computing (3), 1261–1287 (2010).[12] Fidkowski, K. and Darmofal, D. Aerodynamic design optimization on unstructured meshes using the Navier-Stokes equations. AIAA Journal (4), 673–694 (2011).[13] Peter, J., Nguyen-Dinh, M., and Trontin, P. Goal-oriented mesh adaptation using total derivative of aerodynamic functions with respect tomesh coordinates – with application to Euler flows. Computers and Fluids , 194–214 (2012).[14] Belme, A., Alauzet, F., and Dervieux, A. An a priori anisotropic goal-oriented error estimate for viscous compressible flow and applicationto mesh adaptation. Journal of Computational Physics , 1051–1088 (2019).[15] Lions, J. L.
Contrôle optimal de systèmes gouvernés par des équations aux derivées partielles . Etudes mathématiques. Paris: Dunod,Gauthier-Villars, (1968).[16] Sartor, F., Mettot, C., and Sipp, D. Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile.
AIAA Journal (7), 1980–1993 (2015).[17] Morris, M., Mitchell, T., and Ylvisaker, D. Bayesian design and analysis of computer experiments: Use of derivatives in surface prediction. Technometrics (3), 243–255 (1993).[18] Luchini, P. and Bottaro, A. Adjoint equations in stability analysis. Annu. Rev. Fluid Mech. (1), 493–517 (2014).[19] Talagrand, O. and Courtier, P. Variational assimilation of meteorological observations with the adjoint vorticity equation. I: Theory. Q. J. R.Meteorol. Soc. (478), 1311–1328 (1987).[20] Anderson, W. and Bonhaus, D. Airfoil design optimization on unstructured grids for turbulent flows.
AIAA Journal (2), 185–191 (1999).[21] Nielsen, E. and Anderson, W. Aerodynamic design optimization on unstructured meshes using the Navier-Stokes equations. AIAA Journal (11), 185–191 (1999).[22] Nemec, N. and Zingg, D. Newton–Krylov algorithm for aerodynamic design using the navier–stokes equations. AIAA Journal (6), 1146–1154 (2002).[23] Xu. S., Radford, D., M., M., and J.-D., M. Stabilisation of discrete steady adjoint solvers. Journal of Computational Physics , 175–195(2015).[24] Giles, M., Duta, M., and Müller, J.-D. Adjoint code developments using the exact discrte approach. In
AIAA Paper Series, Paper 2001-2596 .(2001).[25] Lu, J. and Darmofal, D. Adaptative precision methodology for flow optimisation via discretization and iteration error control. In
AIAA PaperSeries, Paper 2004-1996 . (2004).[26] Hartmann, R. Adjoint consistency analysis of discontinuous Galerkin discretizations.
SIAM J. Numer. Anal. (6), 2671–2696 (2007).[27] Shi, L. and Wang, Z. Adjoint-based error estimation and mesh adaptation for the correction procedure via reconstruction method. Journal ofComputational Physics , 261 – 284 (2015).[28] Duivesteijn, G., Bijl, H., Koren, B., and van Brummelen, E. On the adjoint solution of the quasi-1D Euler equations: the effect of boundaryconditions and the numerical flux function.
International Journal for Numerical Methods in Fluids , 987–993 (2005).[29] Lozano, C. A note on the dual consistency of the discrete adjoint quasi-one dimensional Euler equations with cell-centred and cell-vertexcentral discretization. Computers and Fluids , 51–60 (2016).[30] Liu, Z. and Sandu, A. On the properties of the discrete adjoints of numerical methods for the advection equation.
International Journal forNumerical Methods in Fluids , 769–803 (2008).[31] Hicken, J. and Zingg, D. Dual consistency and functional accuracy: a finite-difference perspective. Journal of Computational Physics ,161–182 (2014).[32] Stück, A. An adjoint view on flux consistency and strong wall boundary conditions to the Navier-Stokes equations.
Journal of ComputationalPhysics , 247–264 (2015).[33] Stück, A. Dual-consistency study for Green-Gauss gradient schemes in an unstructured Navier-Stokes method.
Journal of ComputationalPhysics , 530–549 (2017).[34] Majda, A. The stability of multidimensional shock fronts. In
Memoirs of the AMS, 275 . Amer. Math. Soc., (1983).
35] Ulbrich, S. A sensitivity and adjoint calculus for discontinuous solutions of hyperbolic conservation laws with source terms.
SIAM J. ControlOptim. (3), 740–797 (2002).[36] Giles, M. B. and Pierce, N. A. Analytic adjoint solutions for the quasi-one-dimensional euler equations. J. Fluid Mech. , 327–345 (2001).[37] Baeza, A., Castro, C., Palacios, F., and Zuazua, E. 2d Euler shape design on non-regular flows using adjoint Rankine-Hugoniot relations.
AIAA Journal (3), 552–562 (2009).[38] Lozano, C. Singular and discontinuous solutions of the adjoint Euler equations. AIAA Journal (11), 4437–4451 (2018).[39] Lozano, C. On the properties of the solutions of the 2D adjoint Euler equations. In Proceedings of EUROGEN 2017, Madrid , (2017).[40] Lozano, C. Watch your adjoints! lack of mesh convergence in inviscid adjoint solutions.
AIAA Journal (11), 4437–4451 (2018).[41] Giles, M. and Pierce, N. Improved lift and drag estimates using adjoint Euler equations. In AIAA Paper Series, Paper 99-3293 . (1999).[42] Todarello, G., Vonck, F., Bourasseau, S., Peter, J., and Désidéri, J.-A. Finite-volume goal-oriented mesh-adaptation using functional derivativewith respect to nodal coordinates.
Journal of Computational Physics , 799–819 (2016).[43] Jameson, A., Schmidt, W., and Turkel, E. Numerical solutions of the Euler equations by finite volume methods using Runge-Kutta time-stepping schemes. In
AIAA Paper Series, Paper 1981-1259 . (1981).[44] Giles, M. and Pierce, N. An introduction to the adjoint approach to design.
Flow, Turbulence, Combustion , 393–415 (2000).[45] Anderson, W. and Venkatakrishnan, V. Aerodynamic design optimization on unstructured grids with a continuous adjoint formulation. Computers and Fluids , 443–480 (1999).[46] Hiernaux, S. and Essers, J.-A. An optimal control theory based algorithm to solve 2D aerodynamic shape optimisation problems for inviscidand viscous flows. In Proceedings of the RTO-AVT Symposium on Aerodynamic Design and Optimisation of Flight Vehicles , (1999).[47] Hiernaux, S. and Hessers, J.-A. Aerodynamic optimization using Navier-Stokes equations and optimal control theory. In
AIAA Paper Series,Paper 99-3297 . (1999).[48] Peter, J. and Dwight, R. Numerical sensitivity analysis for aerodynamic optimization: a survey of approaches.
Computers and Fluids ,373–391 (2010).[49] Hirsch, C. Numerical Computation of Internal and External Flows: The Fundamentals of Computational Fluid Dynamics (second edition) .Butterworth – Heineman. Elsevier, (2007).[50] Goodman, J. and Majda, A. The validity of the modified equation for nonlinear shock waves.
J. Comput. Phys. (3), 336 – 348 (1985).[51] Giles, M. and Pierce, N. An introduction to the adjoint approach to design. In Proceedings of ERCOFTAC Workshop on Adjoint Methods ,(1999).[52] Nguyen-Dinh, M.
Qualification des simulations numériques par adaptation anisotropique de maillages . PhD thesis, Université de Nice-Sophia Antipolis, March (2014).[53] Dwight, R. and Brezillon, J. Effect of approximations of the discrete adjoint on gradient-based optimization.
AIAA Journal (12), 3022–3031(2006).[54] Vassberg, J. and Jameson, A. In pursuit of grid convergence for two-dimensional Euler solutions. Journal of Aircraft (4), 1152–1166(2010).[55] Cambier, L., Heib, S., and Plot, S. The elsA CFD software: input from research and feedback from industry. Mechanics & Industry (3),159–174 (2013).[56] Destarac, D. Far-field / near-field drag balance and applications of drag extraction in cfd, February (2003).[57] Peter, J., Renac, F., Dumont, A., and Méheut, M. Discrete adjoint method for shape optimization and mesh adaptation in the elsA code. statusand challenges. In Proceedings of 50th 3AF Symposium on Applied Aerodynamics, Toulouse , (2015).[58] Anderson, J.
Modern Compressible Flow (third edition) . McGraw-Hill, (2003).[59] Peter, J., Labbé, C., and Renac, F. Analysis of discrete adjoint fields for 2d Euler flows. In
Proceedings of eurogen 2019 . ECCOMAS, (2019).
Appendix A. Additional figures igure A.19: ( M ∞ = . α = o ) Close view of the two 2049 × M ∞ = . α = o , 2049 × δ CL p . δ CL p , δ CL p (equations (13) and equation (14) with ε =1).Figure A.21: ( M ∞ = . α = o , 2049 × δ CL p , δ CL p , δ CL p , δ CL p (equations (13) and equation (14) with ε =1). igure A.22: ( M ∞ = . α = o , 2049 × δ CD p , δ CD p , δ CD p , δ CD p (equations (13) and (14) with ε =1).Figure A.23: (4097 × δ CL p (cid:39) d γ based of three last rows of cells adjacent to the wall. Left : M ∞ = . α = o , right: M ∞ = . α = o .Figure A.24: (2049 × δ R . Left : M ∞ = . α = o , source in (-0.612,-0.127). Right : M ∞ = . α = o , source in (-0.623,-0.020) igure A.25: ( M ∞ = . α = o , 2049 × δ R source term at point (-0.62,-0.02). Discretization of left- and right-hand-side of perturbedCrocco equation (46) ( z component). igure A.26: ( M ∞ = . α = o , 2049 ×2049 mesh) Left: lift adjoint, right: drag adjoint.
Appendix A. Additional figures igure A.19: ( M ∞ = . α = o ) Close view of the two 2049 × M ∞ = . α = o , 2049 × δ CL p . δ CL p , δ CL p (equations (13) and equation (14) with ε =1).Figure A.21: ( M ∞ = . α = o , 2049 × δ CL p , δ CL p , δ CL p , δ CL p (equations (13) and equation (14) with ε =1). igure A.22: ( M ∞ = . α = o , 2049 × δ CD p , δ CD p , δ CD p , δ CD p (equations (13) and (14) with ε =1).Figure A.23: (4097 × δ CL p (cid:39) d γ based of three last rows of cells adjacent to the wall. Left : M ∞ = . α = o , right: M ∞ = . α = o .Figure A.24: (2049 × δ R . Left : M ∞ = . α = o , source in (-0.612,-0.127). Right : M ∞ = . α = o , source in (-0.623,-0.020) igure A.25: ( M ∞ = . α = o , 2049 × δ R source term at point (-0.62,-0.02). Discretization of left- and right-hand-side of perturbedCrocco equation (46) ( z component). igure A.26: ( M ∞ = . α = o , 2049 ×2049 mesh) Left: lift adjoint, right: drag adjoint. igure A.27: ( M ∞ = . α = o , 2049 × igure A.28: ( M ∞ = . α = o , 2049 ×2049 mesh) Left: lift adjoint, right: drag adjoint.