For computation of nonlinear aeroelastic problems, an efficient error estimation and grid adaptation algorithm is highly desirable, but traditional error estimation or grid adaptation do not suffice, since they are insufficiently related to relevant engineering variables and are incapable of significantly reducing the computing time. The dual formulation however, can be used as an a-posteriori error estimation in the quantity of interest. However, derivation of the dual problem, especially the accompanying boundary conditions, is not a trivial task.
This document compares a discrete and analytical ad joint equation method with respect to boundary-condition treatments applied on the quasi-1D Euler equations. Flux evaluation of the primal problem is do ne by a Linearised Godunov scheme. For our future goal, solving ftuid-structure problems, the discrete approach seems preferable

Bijl, H.

Brummelen, E.H. (Harald) van

Duivesteijn, G.F.

Koren, B. (Barry)

For computation of nonlinear aeroelastic problems, an efficient error estimation and grid adaptation algorithm is highly desirable, but traditional error estimation or grid adaptation do not suffice, since they are insufficiently related to relevant engineering variables and are incapable of significantly reducing the computing time. The dual formulation however, can be used as an a-posteriori error estimation in the quantity of interest. However, derivation of the dual problem, especially the accompanying boundary conditions, is not a trivial task. This document compares a discrete and analytical ad joint equation method with respect to boundary-condition treatments applied on the quasi-1D Euler equations. Flux evaluation of the primal problem is do ne by a Linearised Godunov scheme. For our future goal, solving ftuid-structure problems, the discrete approach seems preferable

Koren, B.

Brummelen, van, E.H.

Pure OAI Repository

Comparison of two adjoint equation approaches with respect to boundary-condition treatment for the quasi-1D Euler equations

Duivesteijn, GF

Bijl, H

Koren, B Barry

Brummelen, EH Harald  van

Repository TU/e

CWI's Institutional Repository

International Conference on Adaptive Modeling and Simulation ADMOS 2003 N.-E. Wiberg and P. Diez (Eds) @CIMNE, Barcelona, 2003 COMPARISON OF TWO ADJOINT EQUATION APPROACHES WITH RESPECT TO BOUNDARY-CONDITION TREATMENTS FOR THE QUASI-ID EULER EQUATIONS G.F. Duivesteijn*, H. Bijl, B. Koren, E.H. van Brummelen Faculty of Aerospace Engineering Delft University of Technology P.O. Box 5058, 2600 GB Delft, The Netherlands *email: g.f.duivesteijn@lr.tudelft.nl Key words: adjoint, boundary conditions, quasi-lD Euler, Linearised Godunov Abstract. For computation of nonlinear aeroelastic problems, an efficient error estima-tion and grid adaptation algorithm is highly desirable, but traditional error estimation or grid adaptation do not suffice, since they are insufficiently related to relevant engineer-ing variables and are incapable of significantly reducing the computing time. The dual formulation however, can be used as an a-posteriori error estimation in the quantity of interest. However, derivation of the dual problem, especially the accompanying boundary conditions, is not a trivial task. This document compares a discrete and analytical adjoint equation method with respect to boundary-condition treatments applied on the quasi-1D Euler equations. Flux evaluation of the primal problem is done by a Linearised Godunov scheme. For our future goal, solving fluid-structure problems, the discrete approach seems preferable. 1 G.F. Duive;teijn, H. Bijl, B. Koren, E.H. van Brummelen 1 INTRODUCTION With the increasing reliance on complex nonlinear aeroelastic calculations for complete aircraft, an accurate error estimator for engineering output quantities is of utmost impor-tance. As these computations involve the interaction between flow and structure. these computations may be extremely computing intensive. Use of an efficient grid-adaptation and error-correction technique is required. However, for these applications, traditional error-estimation and grid-adaptation techniques do not suffice, since they are insufficiently related to relevant engineering variables and are incapable of significantly reducing the computing time. A well known strategy for minimising computing costs while achieving a given level of accuracy is grid adaptation. This is done by locally refining the spatial grid in regions where the accuracy of the solution requires this, and by coarsening the grid where the accuracy allows so. A common strategy is to adapt to physical features in the flow, like shocks and boundary layers. However, local refinement of flow features does not guaran-tee a reliable solution for engineering output quantities like lift and drag for aerospace applications. For complex aeroelastic computations, these limitations in a-posteriori error estimation and grid adaptation are even more apparent. Recently, Giles and Pierce [2][3], and Darmofal and Vendetti [61[7] have developed a new method for error estimation and grid adaptation, known as the adjoint-equation approach, shortly: dual approach. The dual formulation is used for the a-posteriori error estimation in the quantity of interest. So far, dual approaches have only been applied to steady flow computations. A promis-ing application area are static and dynamic, non-linear aeroelastic computations. This paper is a first step into that direction. 2 EQUATIONS To illustrate some different ways to implement and use the adjoint method, a well-known test case has been chosen. This test case consists of a convergent-divergent channel, modelled with the quasi-1D inviscid Euler equations. The system of equations is given by: dAF(q) = dA J( ) dx dx q, with With the total energy described as 1 p l 2 E=---+-u "Y - lp 2 ' 2 (1) J(q)- ( D (2) (3) G.F. Duivesteijn, H. Bijl, B. Koren, E.H. van Brummelen and the height A of the channel defined as A(x) -{ ! + sin2 (,rx) (4) the system of equations is balanced. 3 DISCRETE APPROACH The system of equations will be solved by using unsteady finite volume formulation, written as fo a~qdn + in AF(q)dan- fo !: J(q)dn = o. Discretisation of (5) first-order accurate in time yields with the discrete residual operator written as Numerical fluxes are evaluated by using the Linearised Godunov scheme. 3.1 Linearised Godunov scheme (5) (6) (7) To determine the flux vectors across the cell faces, a flux-difference splitting scheme will be used. Besides good capturing of the contact discontinuity and shock waves, the Riemann-based approach is expected to yield robustness and a good boundary-condition treatment. For this reason, a first-order approximation of the full 1D Riemann problem is used to evaluate the flux across a cell face. The flux in or out of a cell is defined as the difference of the individual fluxes at the left cell face and the right cell face. A full derivation of the exact solution of the one- and two-fluid Riemann problem is given in [1]. Regard a cell face ani+1 separating two different states qi and qi+i in ni and Oi+1, as seen in figure 1. In this do~ument, qi and qi+1 will be named q1 and q4 respectively. For our first order flux difference splitting method, the states at nodes i and i + 1 are directly copied on the cell face ani+1, without any interpolation techniques. Recall the Riemann 2 invariants J± in differential form: dJ+ dp along r+, (8) - du+- =0 pa dJ- dp along r-. (9) - du- - =0 pa 3 G.F. Duivesteijn, H. Bijl, B. Koren, E.H. van Brummelen Figure 1: Situation at cell face ani+½. t a) X b) X c) u Figure 2: a) Characteristic curves in physical plane, b) Shock, contact discontinuity and expansion wave in physical plane, c) Hugoniot curve and Poisson curve intersecting at exact solution. The tangents of the two curves intersect at first order accurate solution. In non-homentropic flow without shock discontinuities, the entropy s remains constant when following a particle. However, s may differ for different fluid particles. Due to this non-uniform entropy the Riemann invariants cannot be integrated directly. Using the property that the Riemann invariants J+ and J- stay constant along the characteristic r+ and r- respectively, and linearising (8) and (9) around state 1 and state 4 gives Pl -p2 along rt, (10) U2 -U1 = P1a1 U3-U4 p4 -p3 along f 4. (11) = p4a4 Note that equations (10) and (11) can be interpreted as the tangent of the Hugoniot - and/or Poisson curve in the (u,p)-plane at states 1 and 4. Adding and subtracting equations (10) and (11) yields P2,3 = Pl a1 U1 + p4a4U4 - (p4 - PI) p1a1 + p4a4 p4a4P1 + P1 a1p4 + P1 a1p4a4 ( u1 - u4) (12) (13) with u2 = U3 = u2,3 and P2 = p3 = P2,3· The given solutions u2,3 and P2,3 can be found at the intersection of the tangents of the exact solution curves. This is shown in figure 2c. 4 G.F. Duivesteijn, H. Bijl, B. Koren, E.H. van Brummelen The flux vector fi+1yields 2 (14) The density Pi+l is given by Poisson's relation p = p(p). Hence Pi+½ = P(P2,3) = P2,3· The flux for given c~ntrol volume ni is given by the difference of the two flux vectors at the cell faces: { f:)f:)F dx = Fi+1 - Fi_1. ki X 2 2 (15) Similar expressions can be found when linearising the Osher scheme. For a detailed description of the Linearised Godunov scheme, the reader is invited to read [5]. 3.2 Source term In contrast to the general Euler equations, the momentum equation of the quasi-1D Euler equations contains a source term, in which pressure acts as "propulsion" for the channel, given by the gradient of the geometry. When the pressure of this extra term is evaluated by taking the value at the centre of a given finite volume (implicit or explicit), the discrete solution of the flow field contains an error in comparison with the exact solution. Implicit treatment of the source term in the discrete momentum equation: (Apu):+ 1 = (Apu): + ~ ( Ai-½ (pu2 + p )~-½ -Ai+½ (pu2 + p) ~½ ) + dt ( !: ) i p;+i (16) gives a slight improvement of the computational stability. Moreover, when the pressure in the source term is not evaluated in the cell centre, but over the cell faces instead: (Apu):+1 = (Apu): + ~ ( Ai-½ (pu2 + Pt_½ - Ai+½ (pit2 + p)~+½) + dt ( !: ) i ~ (Pi+½ +Pi-½) k+l , the solution also approximates the exact solution more accurately. 4 DISCRETE ADJOINT METHOD (17) Our goal is to compute an integral quantity, like lift or drag in an efficient way. Suppose the output functional of interest is I= In p(q)dO. 5 (18) G.F. Duivesteijn, H. Bijl, B. Koren, E.H. van Brummelen We use the adjoint-equation method to compute the linear sensitivity of the objective function. For instance, given in [2], when optimising for a single design variable a or change in geometry, the change of the residual function due to a is (19) The influence of the change of a on the discrete output functional can be written as oh ohoqh oh (20) = --+-8a 0% 8a 8a - _ oh ( aRh )-1 oRh + oh (21) 0% OQh Oa Oa TaRh oh Vh Oa + oa' (22) where vh is the discrete adjoint solution of the equation (~:) T Vh = -(~!:) T (23) Another application of the adjoint method, given in [6], is dealing with an accurate solu-tion of the output functional, a solution which is too expensive to be solved on a uniformly fine mesh. The adjoint solution can help to determine where local grid refinement is re-quired to meet the desired accuracy for the computed integral quantity. With Taylor expansion around the coarse grid, we can write for the output functional h(qh) = IH(qH) + [~h] (qh - qf;) + O(h2). qh qf! (24) The solution vector qf! contains an approximated fine mesh solution, which is constructed from a coarse mesh solution. In general, the construction of the fine mesh approximation will be done by high-order interpolation. The fine mesh residual operator Rh(qh) has to be linearised around the coarse mesh: Rh(qh) = Rh(qf!) + [00Rh] (qh - qf;) + O(h2). Qh qf! (25) The Jacobian [~] contains the derivatives of the fine mesh system of equations, % qf! evaluated with the coarse mesh solution. The approximated error vector can be found by isolating (qh - qf!) of equation(25): (26) 6 G.F. Duivesteijn, H. Bijl, B. Koren, E.H. van Brummelen Substituting (26) into (24) yields h(qh) = IH(qH) - ([vh]q,{1) T Rh(qf;), where [vh]q,{1 is the discrete adjoint solution of the adjoint equation 4.1 Adjoint boundary conditions (27) (28) Derivation of adjoint equations from within discrete primal equations has the advantage of not having to worry too much about adjoint boundary conditions. The Jacobian of the residual operator contains the influence of the primal boundary conditions. When taking the transposed Jacobian for computation of the adjoint solution, the adjoint boundary conditions are already included in the system of equations. This can be illustrated in the following way. Writing out the residual operator (7) for volume O; reads where q;_1 , q; and q;+ 1 are the local solution perturbations, and where r; is the residual due to linearisation. This equation gives 3 sub-Jacobians around the main diagonal of %1:, denoted by R,i-l, R;,;, and Ri,i+l· At the boundaries, only 2 sub-Jacobians exist, in which R1,1 and RN,N contain contributions due to implied boundary conditions, e.g. for the outflow at the right hand side of the domain, for the flux for the right cell face FR= FN+l (qN + b.c.). 2 The complete Jacobian looks like: R1,1 R1,2 R2,1 R2,2 R2,3 (31) RN,N-2 RN-1,N-1 RN-1,N RN,N-1 RN,N 7 G.F. Duivesteijn, H. Bijl, 13. Kon'n, E.H. nm Brmnmel,m Tc) evaluate the adjoint solution in a fully numerical way, the Jacobian matrix will be evaluated with a technique, called Automatic Differentiation or Algorithmic ferentiation [4]. A common way of obtaining approximate numerical derivatives of a given function is the divided difference approach, but the main disadvantage of obtaining derivatives in this way is that the method suffers from truncation errors and is prone to reduce the number of significant digits by ½ or ~. In contrast, AD returns derivatives with working accuracy. 5 CONTINUOUS APPROACH Error estimation, based on the adjoint method, can also be performed by the so called continuous method. \Vith this method. described in [3], the analytical adjoint equations will be discretized on the same grid as the primal flow solution. First. ( 1) has to be linearised: I I I dA / · ) I I Lq = AF ( q) · q.r - -d J ( q, · q = r . X (32) The change in the output functional due to small perturbations in the flow solution can be written as I' = { p' · q' dx. Jn (33) The influence of the change in solution on the functional can be determined by the adjoint equation, which can be analytically derived by partial integration: J ( , 1 dA '( ) ') J 'd i 1 • AF · <J.r - x J q · q dx = v • r x, [ , ) '] L J , , dA , , i = v·AF(q -q O - V.r·AF(q)·q +v· dxJ ·qc:r, I I L I I T dA I =[v-AF(q)·q]0 - Jq · [AF(q)] ·t'.r+ [dxJ(q)] •u dx. ( T ) So, the corresponding adjoint equations of (32) are T • [ '( )'T [dAJ'( )] I L v = - AF q 1 . l'.r - dx q . v = P . 5.1 Adjoint boundary conditions (34a) (34b) (34c) (35) Equation (35) can be written by using Jacobians based on the non-conservative flow variables q = (p, u, p) T, so that the adjoint equation becomes, (36) 8 G.F. Duivt'titdjn, H. Bijl. B. Koren, E.H. van Brununelen For this system of adjoint o.d.e\;, boundary and initial conditions have to be defined i order to be able to solve the system. From the derivation of the adjoint equations, show in (34), it follows that the adjoint variables have to be defined in such a way that for qlf · [AF'(qf · v = 0, (37 the dependency of q has been removed. In the next section, three examples will be give in how to derive adjoint boundary conditions for the continuous method. 5.1.1 Subsonic case For the given quasi-1D Euler problem, the boundary conditions of the primal proble1 are defined by p = P·m, u = !Lin at the inflow and p = Pex at the outflow. The nor conservative state vector is written as Looking at small perturbations in the whole domain, the perturbations have to obey tr. boundary conditions as well. At the boundaries, perturbations in the prescribed state are not permitted, so those perturbations have to be zero. q; = 0 } X q~ = 0 q; = 0 ;r Equation (39) written in matrix notation reads: G 0 ~) (~) 1 = q3 (0 0 1) (~) = q3 = -1, = 1. 0 X = -1, 0 X = l, (39~ (391 (41 (4 (4: This results leaves us with one degree of freedom at the inflow and two degrees of freed.a at the outflow. In other words, we are looking for the null spaces to find the missiJ vectors in order to comply (37). Suppose (at x = -1) Null(w1f · [AF'(q)f ·Vin= 0, (4 9 G.F. Duivesteijn, H. Bijl, B. Koren, E.H. van Brummelen where ( 1 0 0) / 0 Wt= 0 1 0 . q = . Rank of w1 is 2, and the kernel has size 1: Null(w1) - Span{w,,} - c, m .Suppose (at x = 1 ) where W2 = ( 0 0 1) · q' = 0. Rank of w2 is 1 and the kernel has size 2: X = -1, X = l. (44) (45) (46) (47) (48) (49) (50) Another example of finding adjoint boundary conditions is given when the boundary conditions for the primal problem are defined by H = Hin, p = Pt;" at the inflow and p = Pex at the outflow. The enthalpy H can be written as H = E + p_ = _i_p_ + !u2 . (51) p 1 -lp 2 When the boundary terms are rewritten in terms of (38), one gets ')'-1 q1 2 2 - m _J.._!Jl + lq2 _ H } q3 + ½q1 q~ = Pt;n X = -1, q3 = Pex X 1. (52) (53) Looking at small perturbations in the whole domain, the perturbations have to obey the boundary conditions as well. At the boundaries, perturbations in the prescribed states are not permitted, so those perturbations have to be zero. ~ ( ~ - ,q~) + q2q2 = 0 } q3 + ½qfo~ + q1q2q2 = 0 q; = 0 10 X = -1, (54a) X 1. (54b) G.F. Duivesteijn, H. Bijl, B. Koren, E.H. van Brummelen Equation ( 54) written in matrix notation reads: X= -1, X = 1. Suppose (at x = -1) where Suppose (at x = 1 ) where W2 = ( 0 0 1) · q' = 0. Rank of w2 is 1 and the kernel has size 2: Null(w,) = Sp&1{w,., w,,) = <;, m + C3 m . (55) (56) (57) (58) (59) (60) (61) (62) (63) Again, multiplying the resulting null vectors with the Jacobian gives the necessary adjoint boundary conditions, which enables us to solve the adjoint equations. 5.1.2 Supersonic case Suppose the primal boundary conditions at the inflow for the supersonic case are given by: M = Min, (64a) H = Hin, (64b) p = Po1n· (64c) 11 G.F. Duivesteijn, H. Bijl, B. Koren, E.H. van Brummelen No boundary conditions are allowed at the outflow. Writing (64) in terms of {38) reads </2 {Gi) 'Y q3 1 2 --- +-qz ') - 1 Q1 2 Rewriting (65) for small perturbations reads l --r:::;r Y'ii </2 0 q3 = Afin, = Hin, Pin· (65a) (65b) (65c) (66a) (66b) (66c) (67} Rank of the matrix is 3 and the kernel 0. This means that for the adjoint equations, no boundary conditions are allowed on the inflow. At the outflow, the kernel contains the identity matrix and gives us 3 boundary conditions for the adjoint equations. Now, the adjoint equation is well posed and can be solved. Numerical results are given in figure 6. The primal problem has been computed, using Linearised Godunov as a flux evaluator. The adjoint solutions have been obtained with the continuous approach. 6 CONCLUSIONS In this document we have compared two adjoint methods. In contrast with the dis-crete method, for the continuous method it is necessary to formulate adjoint boundary conditions. For the discrete method, the adjoint boundary conditions are given implic-itly. Our next step will be the investigation of adjoint based grid adaptation. Some early results have been obtained and look promising. Our goal for the near future is to adapt the spatial grid for some arbitrary output functional, and for a fluid-structure interaction problem. 12 G.F. Duivesteijn, H. Bijl, B. Koren, E.H. van Brummelen ... .------.-----:,:::----..---.=-, -_--, vt .. ·-0 ·---·----····-···-·· ..................... __________________ _ ,.. ,------.-----..----~--:::.,-.-:-.. ,_ Y2-Y3 ........ . ... Figure 3: Adjoint solution for the subsonic (Pt,.,, H;n, Pout) and supersonic (p;n, H;n, M;n) case with h= 3\• REFERENCES [1] G.F. DUIVESTEIJN, Improved capturing of contact discontinuities for two-fluid flows, Note MAS-N0202, CWI, Amsterdam (2002) [2] GILES, M.B. AND PIERCE, N .A., "Adjoint equations in CFD: duality, boundary conditions and solution behaviour", AIAA Paper, 97-1850, (1997) [3] GILES, M.B. AND PIERCE, N.A., "Analytic adjoint solutions for the quasi-one-dimensional Euler equations", J. Fluid Mech., 426, 327-345 (2001) [4] GRIEWANK, A., Evaluating Derivatives, Principles and Techniques of Algorithmic Differentiation, SIAM, Philadelphia (2000) [5] B. KOREN, M.R. LEWIS, E.H. VAN BRUMMELEN AND B. VAN LEER, Riemann-problem and level-set approaches for two-fluid flow computations I. Linearized Go-dunov scheme, Report MAS-R0112, CWI, Amsterdam (2001) [6] VENDETTI D.A., AND DARMOFAL, D.L., "Adjoint error estimation and grid adap-tation for functional outputs: application to quasi-one dimensional flow", J. Comput. Phys, 164, 204-227 (2000) [7] VENDETTI D.A. AND DARMOFAL D.L., "Grid adaptation for functional output: application to two-dimensional inviscid flows", J. Comput. Phys, 176, 40-69 (2002) 13 