Usually we perform numerical simulations in order to control a flow or to improve the design of an object in a flow. Adjoint methods provide a cheap way to calculate gradient information, which can then be used in optimization algorithms. This is an intelligent way to combine CFD with control and design.
Background
Computational Fluid Dynamics (CFD) simulations are often used to improve the design of an object, such as an aeroplane wing, a car, or an inkjet print head. One often sees designers treating CFD like a physical experiment: take the object, calculate the flow, slightly change the object, recalculate the flow, and repeat to an optimal design. This misses a trick.
In the process described above, a new CFD calculation reveals how a single design change influences every aspect of the flow. With a bit of thought, however, the CFD equations can be re-formulated such that a new CFD calculation reveals how every design change influences a single aspect of the flow.
This requires the user to define, ahead of time, what they are trying to achieve. It may be low drag on a wing, or stable flow around a strut. An optimization problem is then formed with this as the cost function and the CFD equations as the constraints. The resultant expression is re-arranged to form the adjoint equations, which show how every parameter or design change will influence the cost function. This provides cheap gradient information, which can then be used in optimization algorithms in order to converge to an optimal design in a handful of iterations. With an appropriate cost function, one can also extract useful physical information from the adjoint solutions.
This describes the basic principles of adjoint-based optimization in Computational Fluid Dynamics Credit: Matthew Juniper, Tejal Shanbhag
Adjoint methods in flow instability
Our first application of adjoint methods was to low density jets and flames. Ultimately, we wanted to calculate the sensitivity of the natural oscillation modes to changes in the flame shape and inlet conditions. In order to do this, we derived the adjoint of Joe Nichols' finite difference low Mach number solver:
This paper describes a derivation of the adjoint low Mach number equations and their implementation and validation within a global mode solver. The advantage of using the low Mach number equations and their adjoints is that they are appropriate for flows with variable density, such as flames, but do not require resolution of acoustic waves. Two versions of the adjoint are implemented and assessed: a discrete-adjoint and a continuous-adjoint. The most unstable global mode calculated with the discrete-adjoint has exactly the same eigenvalue as the corresponding direct global mode but contains numerical artifacts near the inlet. The most unstable global mode calculated with the continuous-adjoint has no numerical artifacts but a slightly different eigenvalue. The eigenvalues converge, however, as the timestep reduces. Apart from the numerical artifacts, the mode shapes are very similar, which supports the expectation that they are otherwise equivalent. The continuous-adjoint requires less resolution and usually converges more quickly than the discrete-adjoint but is more challenging to implement. Finally, the direct and adjoint global modes are combined in order to calculate the wavemaker region of a low density jet.
The low Mach number solver was a good choice because the density is a function of temperature and composition, but not of velocity. In other words, the solver permits hydrodynamic perturbations, but not acoustic perturbations. This means that a pressure correction algorithm can be used with relatively large timesteps and straight-forward boundary conditions.
This was the first time that we had come across the dilemma of whether (i) to write the adjoint of the continuous equations, and then discretize ('continuous adjoint') or (ii) to discretize the continuous equations, and then write the adjoint of the discretized system ('discrete adjoint').
The discretized adjoint equations can be formed in two ways: write the continuous adjoint equations and then discretize (known as the continuous adjoint ) or discretize the direct equations and then take the adjoint of the discretized system (known as the discrete adjoint ). In the DA approach, the direct and adjoint discretized equations have the same truncation errors so rigorous debugging tests can be devised (e.g. by the bug-free code calculating zero to machine precision), which makes it easier to implement. Credit: Matthew Juniper
In the end we did both. Not having anything against which to compare our results, our most important criterion was the ability to debug the code. We wrote the direct solver as a sequence of matrix-vector operations and then coded and debugged the adjoints of those operations one by one. We then debugged the continuous adjoint against the results of the discrete adjoint.
This figure shows the adjoint fields for axial momentum injection, radial momentum injection, and heat injection, for a low density jet (radius 0.5) discharging (at x = 0) into a high density fluid. The cost function is the most unstable eigenvalue. These fields show where the most unstable eigenvalue is most sensitive to the three different types of forcing. This was calculated with the continuous adjoint approach. Note that the fields are smooth. Credit: Gary Chandler Jump to publication (will be at top of next screen)
This figure shows the adjoint fields for axial momentum injection, radial momentum injection, and heat injection, for a low density jet (radius 0.5) discharging (at x = 0) into a high density fluid. The cost function is the most unstable eigenvalue. These fields show where the most unstable eigenvalue is most sensitive to the three different types of forcing. This was calculated with the discrete adjoint approach. Note the oscillations at the entrance plane. Credit: Gary Chandler Jump to publication (will be at top of next screen)
This was also the first time that we came across problems caused by adjoint boundary conditions in Finite Difference codes. Today, we usually express our direct equations in the weak form, so that the direct and adjoint solutions live in the same function space, and use Finite Element or Spectral Element codes. This avoids the easily-avoidable problems. The paper below conducts a thorough examination of this point and contains Matlab examples of good and bad boundary conditions.
Gas turbines and rocket engines sometimes suffer from violent oscillations caused by feedback between acoustic waves and flames in the combustion chamber. These are known as thermoacoustic oscillations and they often occur late in the design process. Their elimination usually requires expensive tests and re-design. Full scale tests and laboratory scale experiments show that these oscillations can usually be stabilized by making small changes to the system. The complication is that, while there is often just one unstable natural oscillation (eigenmode), there are very many possible changes to the system. The challenge is to identify the optimal change systematically, cheaply, and accurately. This paper shows how to evaluate the sensitivities of a thermoacoustic eigenmode to all possible system changes with a single calculation by applying adjoint methods to a thermoacoustic Helmholtz solver. These sensitivities are calculated here with finite difference and finite element methods, in the weak form and the strong form, with the discrete adjoint and the continuous adjoint, and with a Newton method applied to a nonlinear eigenvalue problem and an iterative method applied to a linear eigenvalue problem. This is the first detailed comparison of adjoint methods applied to thermoacoustic Helmholtz solvers. Matlab codes are provided for all methods and all figures so that the techniques can be easily applied and tested. This paper explains why the finite difference of the strong form equations with replacement boundary conditions should be avoided and why all of the other methods work well. Of the other methods, the discrete adjoint of the weak form equations is the easiest method to implement; it can use any discretization and the boundary conditions are straightforward. The continuous adjoint is relatively easy to implement but requires careful attention to boundary conditions. The Summation by Parts finite difference of the strong form equations with a Simultaneous Approximation Term for the boundary conditions (SBP--SAT) is more challenging to implement, particularly at high order or on non-uniform grids. Physical interpretation of these results shows that the well-known Rayleigh criterion should be revised for a linear analysis. This criterion states that thermoacoustic oscillations will grow if heat release rate oscillations are sufficiently in phase with pressure oscillations. In fact, the criterion should contain the adjoint pressure rather than the pressure. In self-adjoint systems the two are equivalent. In non-self-adjoint systems, such as all but a special case of thermoacoustic systems, the two are different. Finally, the sensitivities of the growth rate of oscillations to placement of a hot or cold mesh are calculated, simply by multiplying the feedback sensitivities by a number. These sensitivities are compared successfully with experimental results. With the same technique, the influence of the viscous and thermal acoustic boundary layers is found to be negligible, while the influence of a Helmholtz resonator is found, as expected, to be considerable.
Adjoint methods to reveal active physical mechanisms in a flow
Streamlines (black) and swirl velocity (colour) for the flow around a vortex breakdown bubble (at x = 1.6 to 2.2). This flow is steady but turns out to be unstable. Credit: Ubaid Qadri Jump to publication (will be at top of next screen)
We then applied this code to a uniform density swirling flow in order to detemine the causes of instability around a vortex breakdown bubble. The eigenvalue (natural mode) of the linear stability analysis around a base flow was set as the cost function. We then looked at the influence on the eigenvalue of different types of internal feedback. This revealed two feedback mechanisms at play: (i) a Kelvin-Helmhotz mechanism operating in the shear layer just outside the recirculation bubble and (ii) a feedback mechanism involving conservation of angular momentum just upstream of the bubble.
The flow around the vortex breakdown bubble (at x = 2.2 to 3.0) is unstable to a helical mode with a given natural frequency and growthrate (eigenvalue). Using adjoint methods, we calculate how nine different feedback mechanisms affect this eigenvalue. These mechanisms are for feedback from axial momentum (top row), radial momentum (middle row), and azimuthal momentum (bottom row) causing a force in the axial direction (left column), radial direction (middle column), and aziumthal direction (right column). The influence of feedback in the top row and left column, which is consistent with a Kelvin-Helmholtz mechanism, is greatest at the outer upstream edge of the bubble. Feedback in the bottom-right four frames, which is consistent with a mechanism involving conservation of angular momentum, is greatest around the centreline at the upstream edge of the bubble. Credit: Ubaid Qadri Jump to publication (will be at top of next screen)
At higher swirl numbers, another region of absolute instability appeared but the linear frequency was still determined by the upstream region. The main purpose of this study was scientific: to use adjoint methods to determine instability mechanisms in basic flows.
We then used the same techniques to examine instability in a lifted jet diffusion flame. This revealed two modes: a high frequency oscillation originating from the non-combusting jet and a low frequency oscillation originating from the flame. (The latter is that observed in Juniper, Li, Nichols, Proc. Comb. Inst 32 (1) 1191--1198). This also revealed the best way to control these osillations by controlling the flow.
Streamlines (black) and reaction rate (colours) for a lifted round low density jet diffusion flame. The horizontal axis is the centreline of the jet, which has radius 0.5. The flame sits along the stoichiometric contour (bold black line) and has a distinctive triple-flame structure. This flow is steady but unstable. Credit: Ubaid Qadri Jump to publication (will be at top of next screen)
The lifted round jet diffusion flame is unstable to two types of oscillation: one at high natural frequency and one at low natural frequency. The feedback sensitivity (also known as the structural sensitivity) is large in regions where a generic feedback mechanism has a strong influence on the frequency and growth rate (the eigenvalue). This reveals that the high frequency oscillation (top figure) is generated between the inlet region and the base of the flame, while the low frequency oscillation (bottom figure) is generated in the flame plume. Credit: Ubaid Qadri Jump to publication (will be at top of next screen)
The above analyses were on a cylindrical computational domain and could not easily be extended to more complex shapes. We therefore switched to a Finite Element code that could handle more complex geometry and applied the same techniques:
Adjoint methods in flows with less simple geometry
An X-junction is often used in micro-fluidics to combine three flows. Here, we investigated a symmetric X-junction in which the top and bottom flows have identical parabolic velocity profiles, while the left had flow has a different parabolic profile. Above a critical flow rate, the symmetric flow becomes unstable to a steady asymmetric bifurcation and the exit flow attaches to one of the walls. Credit: Outi Tammisola Jump to publication (will be at top of next screen)
In the X-junction flow, the first bifurcation occurs at a critical Reynolds number and breaks the flow's top-bottom symmetry. This bifurcation can be delayed by symmetric suction/blowing at the wall. This image shows the regions of the flow that, when affected by wall suction/blowing, have most influence on the bifurcation. Credit: Outi Tammisola Jump to publication (will be at top of next screen)
With this code, we considered the influence that blowing and suction at the walls would have on the growth rate of instabilities. One key point was that rounding the walls at the corners increased the transitional Reyolds number greatly. We did not show this formally in that paper, however, because we had not yet put in place the machinery required for shape optimization.
Adjoint methods in turbulent flows with complex geometry
Next, we moved to a more complex geometry and a higher Reynolds number: a gas turbine fuel injector at Reynolds number 4,800.
Hydrodynamic oscillations in gas turbine fuel injectors help to mix the fuel and air but can also contribute to thermoacoustic instability. Small changes to some parts of a fuel injector greatly affect the frequency and amplitude of these oscillations. These regions can be identified efficiently with adjoint-based sensitivity analysis. This is a linear technique that identifies the region of the flow that causes the oscillation, the regions of the flow that are most sensitive to external forcing, and the regions of the flow that, when altered, have most influence on the oscillation. In this paper, we extend this to the flow from a gas turbine?s single stream radial swirler, which has been extensively studied experimentally (GT2008-50278).
The swirling annular flow enters the combustion chamber and expands to the chamber walls, forming a conical recirculation zone along the centreline and an annular recirculation zone in the upstream corner. In this study, the steady base flow and the stability analysis are calculated at Re 200-3800 based on the mean flow velocity and inlet diameter. The velocity field is similar to that found from experiments and LES, and the local stability results are close to those at higher Re (GT2012-68253).
All the analyses (experiments, LES, uRANS, local stability, and the global stability in this paper) show that a helical motion develops around the central recirculation zone. This develops into a precessing vortex core. The adjoint-based sensitivity analysis reveals that the frequency and growth rate of the oscillation is dictated by conditions just upstream of the central recirculation zone (the wavemaker region). It also reveals that this oscillation is very receptive to forcing at the sharp edges of the injector. In practical situations, this forcing could arise from an impinging acoustic wave, showing that these edges could be influential in the feedback mechanism that causes thermoacoustic instability.
The analysis also shows how the growth rate and frequency of the oscillation change with either small shape changes of the nozzle, or additional suction or blowing at the walls of the injector. It reveals that the oscillations originate in a very localized region at the entry to the combustion chamber, which lies near the separation point at the outer inlet, and extends to the outlet of the inner pipe. Any scheme designed to control the frequency and amplitude of the oscillation only needs to change the flow in this localized region.
The large-scale coherent motions in a realistic swirl fuel injector geometry are analysed by direct numerical simulations (DNS), proper orthogonal decomposition (POD), and linear global modes. The aim is to identify the origin of instability in this turbulent flow in a complex internal geometry.
The flow field in the nonlinear simulation is highly turbulent, but with a distinguishable coherent structure: the precessing vortex core (a spiraling mode). The most energetic POD mode pair is identified as the precessing vortex core. By analysing the FFT of the time coefficients of the POD modes, we conclude that the first four POD modes contain the coherent fluctuations. The remaining POD modes (incoherent fluctuations) are used to form a turbulent viscosity field, using the Newtonian eddy model.
The turbulence sets in from convective shear layer instabilities even before the nonlinear flow reaches the other end of the domain, indicating that equilibrium solutions of the Navier?Stokes are never observed. Linear global modes are computed around the mean flow from DNS, applying the turbulent viscosity extracted from POD modes. A slightly stable discrete m = 1 eigenmode is found, well separated from the continuous spectrum, in very good agreement with the POD mode shape and frequency. The structural sensitivity of the precessing vortex core is located upstream of the central recirculation zone, identifying it as a spiral vortex breakdown instability in the nozzle. Furthermore, the structural sensitivity indicates that the dominant instability mechanism is the Kelvin-Helmholtz instability at the inflection point forming near vortex breakdown. Adjoint modes are strong in the shear layer along the whole extent of the nozzle, showing that the optimal initial condition for the global mode is localized in the shear layer.
We analyse the qualitative influence of turbulent dissipation in the stability problem (eddy viscosity) on the eigenmodes by comparing them to eigenmodes computed without eddy viscosity. The results show that the eddy viscosity improves the complex frequency and shape of global modes around the fuel injector mean flow, while a qualitative wavemaker position can be obtained with or without turbulent dissipation, in agreement with previous studies.
This study shows how sensitivity analysis can identify which parts of the flow in a complex geometry need to be altered in order to change its hydrodynamic stability characteristics.
We performed Direct Numerical Simulations and extracted mean flow data and the local eddy viscosity. We then performed linear stability analysis around this mean flow. (This is described in the 'Applied Hydodynamic Stability' project.)
As well as identifying the wavemaker region and sensitivity to steady volume forcing, we could also (in the low Reynolds number case) identify the regions of the injector in which steady suction or blowing would have most influence on the growth rate and frequency of hydrodynamic oscillations.
The Turbomeca gas turbine fuel injector has been extensively studied, both numerically and experimentally. Here, we performed Direct Numerical Simulations (DNS) with a spectral element code at Re = 4,800. This image shows a snapshot of the z-velocity in a slice through the centreline. Credit: Outi Tammisola Jump to publication (will be at top of next screen)
The Turbomecal fuel injector is unstable to a helical instability, which developes nonlinearly into a precessing vortex core. This image shows how suction at the boundary would change the growth rate of the unstable mode. This shows that inflow in the inner channel and at the wall near the separation point is stabilizing. (b) Sensitivity to outwards movement of the domain boundary (shape sensitivity). This shows that changes at two of the corners are influential. Credit: Outi Tammisola Jump to publication (will be at top of next screen)
Adjoint methods in shape optimization
Our next step was to examine the influence of geometry changes. This can be done in two ways: (i) by defining the geometry in terms of parameters (such as the positions of Bezier or B-spline points) and calculating the sensitivity to these parameters, known as parameter-based shape optimization or (ii) by calculating the sensitivity to the positions of the surface gridpoints, known as parameter-free shape optimization. Although parameter-free optimization seems conceptually simpler, it turns out to be technically challenging. Most of our current effort is therefore on parameter-based shape optimization. The paper below shows how to stabilize the flow around a cylinder by altering its shape.
The shape gradients of an instability's growth rate and frequency are derived for an unstable mode calculated from a global stability analysis. These are calculated and interpreted physically for the flow around a cylinder at a Reynolds number of 50. This is a well-known canonical flow, which is often used to discover fundamental behaviour in bluff body flows and to test new numerical techniques. This paper shows that shape deformations affect hydrodynamic oscillations mainly through their influence on the steady base-flow, rather than through their influence on the unsteady feedback mechanism. Deformations that strongly affect the base-flow are shown to strongly affect the frequency and growth rate, as expected. In addition, subtle deformations at the rear of the cylinder are shown to exploit small base-flow changes that have a disproportionately large influence on the growth rate. The physical mechanism behind this is shown to be similar to the well-known phenomenon of ?base bleed?. The method presented in this paper is general and versatile. It provides engineers with gradient information in order to optimize designs systematically. In addition, it provides physical insight, which enables intuitive design changes that would be outside the range of an optimization algorithm or existing geometric parametrization.
This is a simple example of adjoint-based shape optimization. The algorithm changes the shape of the cylinder in order to minimize the viscous drag, while also keeping the flow marginally stable. Credit: Jack Brewster
We have applied this to shape optimization of the cyclonic separator in a bag-less vacuum cleaner:
In many engineering applications it is desirable to have a steady flow field that is stable to perturbations. This can be described through a global stability analysis. This is an eigenvalue problem which gives a series of mode shapes (the eigenfunctions) and their growth rates and frequencies (the eigenvalues).
This thesis calculates and interprets the shape sensitivity of the eigenvalue of a global stability analysis. The shape sensitivity allows quick calculation of the change in a mode?s growth rate and frequency when the flow geometry is changed. The shape sensitivity is calculated using an adjoint method. This is derived both for flows governed by the incompressible Navier?Stokes equations and also for flows governed by the URANS equations using the Spalart?Allmaras and k-omega turbulence models. Using the laminar and turbulent flow past a cylinder as a model problem, the shape sensitivities of the laminar and turbulent vortex shedding modes are presented. For both cases, control of the eigenvalue by shape deformation is shown to occur primarily by changing the steady base-flow and not by changing the mode?s unsteady feedback. A technique for finding the deformations with the greatest influence on the base-flow is then demonstrated. For the laminar flow, this shows that the shape sensitivity of the growth rate and frequency is primarily due to a single deformation that causes large widespread base-flow changes. This deformation increases the growth rate but decreases the frequency, leading the shape sensitivities to have similar shapes but opposite signs.
The thesis then applies this framework to the cyclone separator. This is an industrial application which has been shown to have a helical instability that reduces the cyclone's performance. Helical and double-helical modes with similar Strouhal numbers to those seen in experiment are found. As with the cylinder, control of the eigenvalue by shape deformation is shown to occur primarily by changing the steady base-flow and not by changing the mode?s unsteady feedback. The shape sensitivities are shown to be concentrated at the cyclone's cone tip and vortex finder. Small changes of these parts of the geometry are shown to cause large widespread changes to the base-flow. Finally, the shape sensitivities are used in a gradient based method to show that small changes in the cyclone geometry can significantly reduce the growth rate of the unstable modes. This is done in conjunction with a separation metric to show that the stability of the flow can be improved without reducing separation performance.
The techniques demonstrated in this thesis for finding the shape sensitivity of the eigenvalue and for identifying influential deformations can be easily applied to a wide range of different flow geometries, governing equations and cost functions. This highlights the benefits of the use of adjoint-based methods in engineering design.
We have applied shape optimization to inkjet print heads:
In drop-on-demand injket printheads, ink is pumped steadily through small channels, each of which contains an actuator and a nozzle. When an actuator pulses, a droplet is forced through the nozzle, after which acoustic oscillations reverberate within the channel. Manufacturers would like to damp the residual reverberations, without increasing the pressure drop required to drive the steady flow. In this paper we show that this can be achieved by constricting the channel where the acoustic velocity is largest and enlarging the channel where the acoustic velocity is smallest. This increases the viscothermal dissipation of the acoustics without changing the viscous dissipation of the steady flow. We separate the compressible Navier--Stokes equations into equations for a steady flow with no oscillations and equations for oscillations with no steady flow. We define two objective functions: the viscous dissipation of the base flow and the viscous dissipation of the oscillations. We then derive the adjoints for both sets of equations, and obtain expressions for the gradient of each objective function with respect to boundary deformations in Hadamard form. We combine these with a gradient-based optimization algorithm, incorporating constraints such as the shapes of the actuator and nozzle. This algorithm quickly converges to a design that has the same viscous dissipation for the steady flow but a 50% larger decay rate for the oscillating flow. We show that this design is nearly optimal. It is a shape that inkjet manufacturers, using physical insight and trial and error, have not yet considered. We also show how the adjoint fields provide physical insight into the mechanisms affecting each objective function. The main requirements of this method are that the steady flow Mach number and oscillating flow Mach number are small, that the oscillations have relatively small amplitude, and that dissipation is dominated by thermo-viscous mechanisms. These requirements are often satisfied in microfluidics, so the method in this paper could be applied to many other applications.
We have applied shape optimization to stabilize thermo-acoustically unstable combustors:
We construct a low order thermoacoustic network model that contains the most influential physical mechanisms of a thermoacoustic system. We apply it to a laboratory-scale longitudinal combustor that has been found to be thermoacoustically unstable in experiments. We model the flame, which is behind a bluff body, by a geometric level set method. We obtain the thermoacoustic eigenvalues of this configuration and examine a configuration in which six eigenmodes are unstable. We then derive the adjoint equations of this model and use the corresponding adjoint eigenmodes to obtain the sensitivities of the unstable eigenvalues to modifications of the model geometry. These sensitivities contain contributions from changes to the steady base flow and changes to the fluctuating flow. We find that these two contributions have similar magnitudes, showing that both contributions need to be considered. We then wrap these sensitivities within a gradient-based optimization algorithm and stabilize all six eigenvalues by changing the geometry. The required geometry changes are well approximated by the first step in the optimization process, showing that this sensitivity information is useful even before it is embedded within an optimization algorithm. We examine the acoustic energy balance during the optimization process and identify the physical mechanisms through which the algorithm is stabilizing the combustor. The algorithm works by, for each mode, reducing the work done by the flame, while simultaneously increasing the work done by the system on the outlet boundary. We find that only small geometry changes are required in order to stabilize every mode. The network model used in this study deliberately has the same structure as one used in the gas turbine industry in order to ease its implementation in practice.
Ensembles of geophysical models improve prediction accuracy and express uncertainties. We develop a novel data-driven ensembling strategy for combining geophysical models using Bayesian Neural Networks, which infers spatiotemporally varying model weights and bias, while accounting for heteroscedastic uncertainties in the observations. This produces more accurate and uncertainty-aware predictions without sacrificing interpretability. Applied to the prediction of total column ozone from an ensemble of 15 chemistry-climate models, we find that the Bayesian neural network ensemble (BayNNE) outperforms existing methods for ensembling physical models, achieving a 49.4% reduction in RMSE for temporal extrapolation, and a 67.4% reduction in RMSE for polar data voids, compared to a weighted mean. Uncertainty is also well-characterized, with 91.9% of the data points in our extrapolation validation dataset lying within 2 standard deviations and 98.9% within 3 standard deviations.
