Restricted snakes volume of solid (RSVS): A parameterisation method for topology optimisation of external aerodynamics

An approach for aerodynamic shape optimisation is derived which is capable of handling topological design changes as well as detailed surface control. The technique applies a material distribution, or volume of solid approach where design variables specify a volume fraction of solid on a fixed mesh. To convert this data to a solid surface, a contour is constructed around the volumes by moving points on the surface until the final shape satisfies those specified volumes. The objective of this construction procedure is to minimise the surface length, subject to the preset volume constraints. As a result, the method reproduces circular arcs exactly. Shape function analysis is then used to explore the theoretical behaviour of the parameterisation, and to prevent oscillatory surfaces from forming, thereby ensuring good optimiser convergence. The method is extended to allow for anisotropic refinement of the parameter mesh. Final test cases include geometric fitting of arbitrary shapes, as well as drag minimisation of topologies in supersonic flow, and show the parameterisation is able to explore single and multi-body aerodynamic


Introduction
Increases in computational power and improvements in computational fluid dynamics (CFD) tools have created the possibility of using CFD-based optimisation in industrial design. By allowing a systematic and unbiased exploration of a design space, optimisation methods can be used to expand a designer's understanding of the problem being tackled, allowing better overall aerodynamic performance. As designers look to improve performance, aircraft manufacturers are turning increasingly to numerical optimisation. Frameworks for aerodynamic optimisation require the integration of parameterisation methods, mesh generators and flow solvers with optimisation methods. The tendency in this has been to use a modular approach by integrating established modelling and CFD packages with existing optimisers.
The complexity of parameterisation arises from the different origins of optimisation methods and CFD processes. Optimisation methods are mathematical algorithms devised to find the extrema of functions, and have rigorous mathematical underpinnings, while CFD originated from the need to evaluate the aerodynamic properties of potential designs. The translation of the mathematical for-mulations used by optimisers into the geometric designs used by CFD is a complex problem with implications for the efficiency and effectiveness of optimisation frameworks. Parameterisation methods for aerodynamics need to be compact while not artificially limiting the geometric shapes that can be represented [1,2] . This focus led to aerodynamic optimisation methods capable of efficiently handling small surface changes, using 10s to 100s of design variables in 2 dimensions and 100s to 10 0 0s in 3 dimensions. While the compactness of aerodynamic parameterisations improves the convergence of optimisers, it has come at the cost of a more limited capability for creating topological changes.
In structural design the benefits of exploring different topologies are key to generating efficient structures. The field of numerical structural topology optimisation (STO) has been an active field of research for the last 30 years and it has recently seen industrial application on the Boeing CH-47 Chinook and the Airbus A380; it allowed a weight reduction of 17% of underfloor beams compared to a conventional structural optimisation method [3] on the CH-47 and weight reduction of the leading edge droop ribs on the A380. This effort in the finite elements (FE) community has led to parameterisation methods able to represent complex topologies with a homogeneous set of design variables [4][5][6] . Recent progress in STO has culminated in the numerical optimisation of an entire Boeing 777 wing under aerodynamic loads with 1.1 billion degrees of freedom by Aage et al. [7] . The algorithm used in that work was capable of building features 40 mm long in the 27 m half-span, resulting in very detailed internal structures resembling bone and beak structures in nature.
The justification for topological optimisation is straightforward in structural applications, from truss space-frames to honeycomb designs, there are a wide range of possible engineering structures; furthermore a structural member's impact is readily summarised to a set of interactions at its boundary. The possibility to reduce designs to a set of external interactions and the Lagrangian formulation of CSD solvers facilitates the implementation of structural topological optimisation within existing designs.
There is no such separation in aerodynamics; the aerodynamic shape is intrinsically linked to the rest of the design by its need to be supported by an underlying structure. This means that aerodynamic topological optimisation of an entire aircraft or wing is unlikely to be a reality in the near or medium term. However, there is scope for the aerodynamic topological optimisation of local features; topological optimisation of wing tips would allow feathered or split winglets of the type seen on the MD-11 and Boeing 737-MAX to be explored. No current optimisation framework for external aerodynamics supports the exploration of topological changes, because none of the parameterisation methods commonly in use can represent different topologies with a homogeneous set of design variables.
An effective topological aerodynamic optimisation framework offers the possibility of radically new designs. Applications to Formula 1, unmanned aerial vehicles, commercial strut-braced wing design and internal engine design could offer significant improvements in performance. This paper presents the development of a parameterisation method which can handle topology changes while maintaining a compact design space, allowing the exploration of new aerodynamic optimisation problems.

Existing aerodynamic and structural optimisation methods
Earlier developments in the field of parameterisation for aerodynamics have yielded a wealth of different methods for the representation of aerodynamic designs. Parameterisation methods can be separated broadly in two categories: constructive and deformative methods. Constructive methods define completely the geometry from the set of design variables; these include B-Spline and polynomial interpolation [8] in general, and CST [9] and PAR-SEC [10] in particular. Deformative methods by comparison define a set of modifications to a baseline geometry; notable among these are the Hicks-Henne bump functions [11] , Singular Value Decomposition (SVD) deformation modes [12,13] and Free-Form Deformation (FFD) methods [14,15] . While most parameterisations presented here can be extended to three dimensions, their capability varies widely. In three dimensions a common approach is to use FFD deformation methods as these can be adapted to work directly on an existing mesh.
Previous systematic investigations by Vassberg et al. [1,16] have highlighted the impact of dimensionality on the drag minimisation of a standard test case, showing the importance of geometric flexibility while maintaining a compact set of design variables. Work by Castonguay and Nadarajah [17] , and more recently by Masters et al. [18,19] has compared the impact of established parameterisation methods on geometric flexibility, pressure distribution recovery and optimal drag results. These studies show that effective parameterisation methods require few design variables while still maintaining smooth control of the profile. Smooth control is achieved when a small change in the numerical representation leads to a similarly small change in the represented geometry. This requirement results from the expense associated with converging optimisers in large design spaces, balanced against the need to not artificially restrict the scope of geometries that can be represented [1] .
Most aerodynamic parameterisation methods to date have focused on producing smooth designs with small numbers of design variables (in the 10s to low 10 0 0s). One key restriction that affects nearly all established parameterisation methods is the inability to transition between topologies, so split or multi-body configurations cannot be explored. This work presents an aerodynamic parameterisation that does provide this level of topological flexibility.
In structural topology optimisation homogenisation and level set methods have been used to tackle complex problems in two and three dimensions; however these structural methods have limitations in terms of their application to aerodynamics. The first methods developed for STO were homogenization methods ; these rely on the segmentation of the design domain into squares in which the density of a material can be varied to change a design's weight and local load carrying ability. By affecting directly the density of the material in the discretization of the structural solver homogenization methods do away completely with the need for an explicit representation of a profile. These works led to the development of the solid isotropic material with penalisation (SIMP) method [20] , the most widely used STO procedure. Homogenization methods do not maintain a representation of the outer boundary of the shape and instead rely on direct interaction with the structural solver.
The main alternative to homogenization are level-set methods (LSM) introduced by Wang et al. [21] . In these methods the structural profile is represented by the level set of a parametric function. These methods were shown to be very competitive and solve some of the shortcomings of homogenization methods [21] . Level sets methods include a wide range of approaches for the definition of the level set function; each of these choices affects the behaviour of optimisation processes [22] . However these methods have in common the implicit definition of the profile and can rely on three mechanisms for change: boundary profile variations; functional parameter variations; and topological variations. To be effective these methods rely on very close integration with the optimisation method, usually through adjoints. Originally this was done by using the Hamilton-Jacobi updates to propagate the boundary of the profile as a moving front [22] . More recently, mathematical programming approaches such as sequential quadratic programming (SQP) and the method of moving asymptotes (MMA) have become popular as well as global, gradient free approaches [22] . For a complete overview of the field of structural topology optimisation the reader is referred to the comprehensive review by Deaton and Grandhi [4] .
Both homogenization methods [23,24] and level set methods [25] have been adapted to fluid topology optimisation in two and three dimensions. Modelling the Stokes equations and incompressible flows at low Reynolds numbers, solved by finite element elasticity solvers, these methods have yielded good results on the optimization of micro-fluidic devices and channel flows [24,26] . Lattice-Boltzmann methods have been used to tackle some very low Reynolds number problems [27,28] . Most of these methods rely on derivatives with regard to the porosity of the material to drive the evolution of the topology. Recent advances have seen broader ranges of turbulent incompressible flows being tackled [29,30] . These methods do not maintain a smooth and crisp fluid boundary, limiting their use for compressible aerodynamics which use solutions to the Euler and compressible Reynolds averaged Navier-Stokes (RANS) equations at high Reynolds numbers. One notable exception is the cutFEM method by Villanueva et al. [31] , where a LSM is used to parameterise shape and topology at very low Reynolds number.  Some progress towards a topologically flexible compressible aerodynamic optimisation has been achieved by Hall et al. [32,33] . This method relies on material distribution, or volume of solid (VOS), to generate the external geometry of an aerodynamic body. An example of this type of design space is shown in Fig. 2 (a). The VOS approach is inspired from volume of fluid methods used in multi-fluid simulations. It is, in optimisation terms, conceptually similar to density approaches: it defines explicitly regions of space which are full or empty based on a predefined grid. In each cell of this grid the fraction of that cell which must be inside the profile is specified by a value between 0 and 1, this allows the parameterisation to be understood intuitively by a designer. The VOS method by Hall et al. [32] uses this information to generate a smooth level set function from which a contour that approximately matches the VOS is extracted. While effective on cases where topological flexibility was required, it required more design variables compared to other aerodynamic parameterisation methods for similar geometric accuracy.

Development of the aerodynamic topology optimisation framework
These observations show that the development of a topological aerodynamic optimisation framework has the potential to deliver further improvements in both conventional and future aircraft configurations. This paper presents development of the r-snake volume of solid (RSVS) method, an aerodynamic parameterisation that supports topological change while still performing efficiently on typical problems. To be useful, the RSVS needs to fit into cur-rent modular aerodynamic frameworks, it must have: a sufficiently compact and smooth design space; be compatible with Eulerian finite volume CFD approaches; support adjoint gradients and be extensible to three dimensional problems. The framework developed in this paper is summarised in Fig. 1 . To ensure compatibility with the optimisation methods already shown to be effective in aerodynamic optimisation the set of design variables needs to be homogeneous; that is to say all design variables must be of the same type. This precludes the use of traditional aerofoil parameterisation methods with additional variables explicitly controlling the topology of the geometry.
The RSVS builds upon the volumetric aerodynamic parameterisation by Hall et al. [32,33] which was an early topologically flexible parameterisation for external aerodynamics. Like the parameterisation of Hall et al., the RSVS uses volume of solid (VOS) design variables to control profile shape and topology, these were kept as they provide intuitive handling of topology change. However an effective contour generation method has been developed to improve the geometric behaviour of the parameterisation. The main challenge in this type of parameterisation is the translation of the design variables into profiles suitable for CFD analysis. This paper details how the RSVS allows to go from a VOS design space ( Fig. 2 (a)) to a specific two body profile built for a set of VOS design variable values ( Fig. 2 (c)).
In volume-based parameterisation the segmentation of volumetric information is done through a Cartesian grid, this means the design variables are best understood by a designer as grey-scale images on an underlying mesh ( Fig. 2 (a)). This observation highlights the similarity between the parameterisation of geometries from volume information and the field of contour extraction in image analysis. Image segmentation, and medical image segmentation in particular, pose many of the same challenges as the volumetric parameterisation method considered earlier. The recovery of complex closed contours of arbitrary topology with limited computational expense is one that has been explored by the medical imaging community for the last 20 years. A class of methods for building such profiles that has seen significant and promising use is that of active contour methods [34,35] . These methods rely on explicit vertex marching until the contour meets internal and external forcing conditions. Restricted snakes (r-snakes) developed by Kobbelt and Bischoff [36] are a type of parametric active contour designed to handle topology changes efficiently. Section 2 shows how r-snakes are used in the RSVS to generate profiles of suitable aerodynamic quality that respect the values of VOS design variables. The shape of the r-snake is driven by a set of equations that were found to have desirable smoothness properties, these are presented in Section 3 . This process was used to generate the profile in Fig. 2 (c). This approach to contour generation and representation has the benefit of being fully explicit, bringing it in line with current established parameterisation methods.
In addition to the r-snake volume of solid parameterisation, significant work has been carried out to understand and improve the performance of the VOS design variables used by the RSVS. Design variable smoothing processes and refinement are developed in this work to improve the behaviour of optimisation frameworks relying on the new parameterisation method. Multi-level approaches to parameterisation by Anderson and Aftosmis [37] and Masters et al. [38] have shown their ability to accelerate and improve the performance of underlying optimisation frameworks. A similar hierarchical method is developed for the RSVS parameterisation in Section 5 . This multi-level approach allows significant performance improvements on the basic RSVS implementation in geometric and aerodynamic optimisations while removing some of the expert knowledge required when setting up new optimisation cases.
To validate these developments geometric and aerodynamic design cases were explored using the RSVS parameterisation method with refinement. Integration of the parameterisation with optimisation methods is presented in Section 6 . These cases aimed to exploit both the topological flexibility of the parameterisation as well as its compact and smooth formulation. In Section 7 , the geometric inverse design of single and multi-body airfoils is tackled; in Sections 8 and 9 results for the drag minimisations of areaconstrained single and multi-body profiles in supersonic flow are shown.

Geometry and topology generation using restricted snakes
The role of the parameterisation method is to provide an interface between an optimisation method and a solver to form a shape optimisation framework. Efficiency and flexibility of shape optimisation frameworks is limited by the geometric capability of the parameterisation method. This section presents how the r-snake volume of solid (RSVS) parameterisation translates sets of volume fraction design variables specified on a fixed grid into closed contours of varying topology. For optimisation frameworks to exploit the RSVS efficiently, this process must reliably produce smooth features at a resolution below the grid on which VOS values are defined.
To achieve this level of smooth control, the RSVS profile is defined as the closed contour of minimum arc-length that will match the volumes of the design variables; it is built using a restricted snake (r-snake). The r-snake (e.g. in Fig. 5 ) is a method for "vertex marching" which allows efficient topology handling and is tolerant of any layout of VOS design variables. The r-snake is a type of parametric active contour originally developed by Kobbelt and Bischoff [36] . This section develops the integration of the rsnake with the RSVS condition of minimising the arc-length under volume constraints. This condition was found to reliably produce smooth profiles enabling a compact parameterisation. Later sections explore the analytical properties of RSVS equations to show that the stated smoothness and compactness targets have been achieved. The RSVS process is summarised in Fig. 3 .

Formulation of the RSVS
One of the difficulties in designing a parameterisation with topological flexibility is to maintain smooth control close to topology changes, as these are geometrically discontinuous regions of the design space. To define a set of VOS variables a grid is superimposed on the design space, where the design variables become the fraction of each cell within a geometry built from this information. This process is shown for a simple grid in Fig. 4 . This parameterisation procedure provides intuitive handling of topology change without maintaining explicit control of it. It is important that topology is not controlled explicitly as this would lead to a severely discontinuous design space which would not be usable with many of the traditional local and global optimisers used for aerodynamic optimisation.
The VOS design variables do not include in themselves rules for building a profile. These rule must generate profiles which are continuous and smooth, allow features smaller than the VOS design variables, and be indifferent to the type of grid they are being applied to. This last requirement opens up the possibility of using non-square grids for improved flexibility and compactness of the method. The rule must also be extensible to the generation of water-tight surfaces in three dimensions with minimal modification.
The condition used to define the RSVS is minimisation of the profile length, with the constraint that the area enclosed by the contour within each design cell must exactly match the value for  of the VOS. The mathematical formulation of this problem is given in Eq. (1) . This system is analogous to the effect of a tensile force "shrink-wrapping" the required VOS in each cell; the benefit is it allows for smooth profiles in most cases but can also recover sharp corners where the VOS requires it. This formulation also has a natural extension in three dimensions as the surface area minimisation under volume constraints of a geometry.
In the expression above y is the closed profile built by the RSVS, m is the number of VOS cells in the design space, a j the value of the VOS and C j the outer boundary of the j{th} cell. These are represented graphically in Fig. 4 (b). The VOS is taken as a constraint on the area enclosed in both the profile and each cell. The next sections detail how this mathematical program can be solved using restricted snakes to produce an effective shape and topology parameterisation method.

Marching of the r-snake
To build the RSVS parameterisation method the r-snake must be evolved until it solves the length minimisation problem specified in Eq. (1) . To allow a high degree of geometric flexibility with few design variables, features need to be recovered below the resolution of the VOS grid. The restricted snake is a vertex marching procedure where the control points (called snaxels) are constrained to move on a predefined grid, as a consequence this snaking grid controls the number of snaxels and the resolution of the contour. By marching the snake on a grid finer than the VOS grid, smooth features below the resolution of the volume design variables can be recovered.
To drive the position of the r-snake the original continuous length minimisation problem ( Eq. (1) ) is discretised in terms of the r-snake and snaxel variables into the mathematical program in Eq. (2) . This discretization process needs six properties from the r-snake geometry and the snaxel positions. The first three of these properties are part of the snaking algorithm; the last three properties of the snaxels are derived from connectivity and grid information, and are needed for the implementation of the discrete length minimization problem. These properties are: the snaxel index ( i ), used to reference it in all operations; the normalised position along an edge ( d i ∈ [0, 1]); the scalar velocity along that edge ( v i ∈ R ); the snaxel position in Cartesian coordinates ( p i ); the direction of travel of the snaxel ( g i ) and the vertex of origin ( g i,1 ); the normal vectors to the preceding and following edges ( n i and n i +1 ). These properties are represented graphically on Fig. 5 In Eq. 2 d is the column vector of all snaxel distances d i , a is the column vector of target volumes a j in each VOS cell and A ( d ) is the current volume in each VOS cell contained by the r-snake. The normalised snaxel positions ( d ) are used as the design variables of the length minimisation problem. This formulation has the benefit of being very general, it can be tackled on an arbitrary volume grid with any underlying snaking grid with any optimisation method. This generality guarantees a high degree of flexibility in the range of shapes that can be represented. The following sub-sections show how this problem can be solved efficiently by using a Newton step sequential quadratic programming (SQP) procedure. The availability of analytical gradients for the profile length condition (the objective function) and the volume information (constraints) means a gradient based method may be applied efficiently.
To solve this surface length minimisation a method was required that would converge in few iterations and function evaluations. For this reason a sensible choice is to use a gradient based method. The availability of analytical first and second derivatives means that sequential quadratic programming (SQP) is a viable option. A damped Newton step defined from a quadratic approximation to the full mathematical program is used to specify the snaxel velocities. The full derivation of the Newton step SQP is originally presented in Boggs and Tolle [39] . Only the final velocity update formula used in the RSVS is presented below ( Eq. (3) ).
The change in distances k +1 d is used as the velocities ( v i ) of the snaxels, letting the snaking process handle damping and connectivity changes. Calculation of the derivatives necessary to the evaluation of Eq. (3) is done analytically. The derivation of these is beyond the scope of the main body of this paper but is useful for the implementation of the RSVS. It is therefore presented in Appendix A . Validation of this process is presented in Section 4.1 .

Topology initialisation and evolution
While the RSVS rules specified in the previous section define desirable properties for the profile and a method to evolve an existing profile, they do not specify a starting geometry. External aerodynamic optimisation is usually concerned with the design of outer surfaces, for this reason the r-snake is initialised at the outer boundary of non-empty VOS cells. An example of this type of initialisation is shown in Fig. 6 , by the outer red profile. This approach to initialisation has the benefit of always being defined and has an intuitive behaviour: it is similar to a force shrink wrapping the VOS design variables.
The benefit of using the restricted snake is that the topology can be modified if the volume fractions require it. This process is again illustrated in Fig. 6 by the intermediate blue profiles. The dashed contour shows the r-snake before topology cutting and the solid blue line after the topology cut.
The topology change is handled efficiently by the restricted snake. The change in the geometry around the cut (shown in the close-up in Fig. 6 ) is due to the r-snake algorithm removing invalid snaxel connections as specified by Bischoff and Kobbelt [40] . To maintain the integrity of the profile the r-snake algorithm limits the possible connections of a snaxel with it's neighbours. The rules as developed by Kobbelt and Bischoff [36] are: no 2 connected snaxels can be on the same edge; snaxels must travel out of the profile. When two snaxels meet the profile connectivity is altered to by-pass them in what results in a change of profile topology. The connectivity rules are then applied removing the invalid connections that were generated. This process is detailed in [36] .

Extension to 3-dimensions
One of the benefits of the RSVS is that all the elements are naturally extensible to three dimensions. In three dimensions, the governing equation becomes the minimisation of surface area under volume constraints of the generated geometry. Parametric active contour with topological flexibility have also been used before for the segmentation of 3-dimensional medical images [34] , and the extension of the restricted snake to a surface is possible. Finally an SQP can still be used to update snaxel positions thanks to the differentiability of the problem in 2 and 3 dimensions. The challenge in 3 dimensions is that the SQP system solved using Eq. (3) is considerably larger (thousands of snaxels) making the use of direct single core solvers slow. Fortunately, the matrices of the system are sparse allowing the efficient use of parallel solvers and iterative methods.

Behaviour of the restricted snake volume of solid (RSVS) parameterisation
Performance of aerodynamic optimisation frameworks is highly dependent on the behaviour and flexibility of the parameterisation method. To understand how the RSVS will perform in this context it is useful to explore the properties of the curves generated by this new method. The discrete constrained length minimisation which governs the behaviour of RSVS profiles can be expressed for an analytical continuous profile. This analytical formulation to the problem can be explored using calculus of variations and this section shows that piecewise continuous circular patches are the solution to the analytical RSVS problem. This result is validated in Section 4.1 by showing the geometric convergence of practical profiles with increasing snaxel densities compared to an analytical solution. Also of interest in shape optimisation is the response of profiles to small changes of the design variables; the analytical formulation of the length minimisation allows the calculation of the shape response due to small changes in the volume constraints which are the design variables for shape optimisations that use RSVS. This section shows that for well parameterised profiles the response is a quadratic spline.
The formulation of the RSVS presented in Eq. (1) can be analysed using calculus of variations to arrive at analytical expressions for the profiles generated. This analysis shows that the curves generated are continuous splines made of one arc of circle in each design cell. These curves can themselves be represented as second order non-uniform rational b-splines (NURBS). To make the derivation of these curves more straightforward it is useful to consider the simplified case of a curve minimising length between two points of coordinates ( b 1 , 0) and ( b 2 , 0) with a single volume constraint of value a . This minimisation problem needs to be solved for the expression of curve y ( x ). The optimum is the solution to the Euler-Lagrange equation of this problem which is expressed in Eq. (4) . Eq. (4) can be solved for y and then for positive values integrated into Eq. (5) .
Constants c 0 and c 1 are integration constants, λ is the Lagrange multiplier, these variables are to be chosen based on limit conditions. The equation for y ( x ) is found to be the equation of a circle in Cartesian coordinates. To understand the effect of the various integration coefficients it is useful to rearrange Eq. (5) into the canonical form of the circular equation.
From Eq. (6) it becomes clear that c 0 and c 1 control the x and y positions, respectively, of the centre of a circle of radius 1/ λ. Exactly the same reasoning that was applied to a single volume constraint can be extended to a multi-constraint problem corresponding to a full RSVS profile. The many constraint problem shows that the RSVS profile is a piecewise continuous curve with each patch an arc of circle. C 0 and C 1 continuity of this piecewise profile come out of the assumption behind the Euler-Lagrange process ( Eq. (4) ) used for solving the RSVS governing equation.
One of the benefits of the length minimisation problem resulting in a set of patched arcs of circles is that it can be represented exactly by NURBS of degree 2. The ability of the RSVS to generate curves of this class is likely to enhance its usability at early design stages, and allows existing work on NURBS to be leveraged in the design of refinement procedures and equivalences with other geometry generation tools. To guarantee smooth pressure distributions, aerodynamic shapes are expected to be curvature (G2) continuous. In optimisation this needs to be traded off with the geometric flexibility and generality of the parameterisation: wings and aerofoils while mostly smooth present discontinuities at the trailing edges that must be represented. The derivative continuity of the RSVS allows the representation of sharp corners and smooth shapes. Continuity of higher order derivatives can be achieved by the optimiser naturally producing smoothly distributed design variables. If this additional levels of smoothness is required smoothness constraints can be applied on the design variables to achieve much the same effect as a smoother parameterisation [41] .
A similar process can be followed to analyse the response of the RSVS to a small change in volume fraction. This is a useful property as the performance of gradient-based aerodynamic optimisation frameworks depends on the geometric change caused by small changes in design variables. While not presented here the derivation shows that, for a profile with low curvature in the parameter cells, the response will tend to a C 1 continuous quadratic spline. Practical uses of RSVS result in profiles with low curvature in each design cell which means that the response of the RSVS to a small disturbance will be close to a quadratic spline, a well understood class of functions.

Parameterisation results
This section presents and discusses profiles and behaviour of the RSVS parameterisation method. The profiles shown in this section were all designed with fixed volume fraction values to test the capability of the parameterisation method; profiles generated by shape optimisation frameworks where VOS values are changed iteratively as design variables are shown in Sections 7 and 8 . Section 4.1 presents validation of the RSVS parameterisation while sub-Section 4.2 shows results for the geometric recovery of common aerofoil sections.

Validation tests
This section presents some of the validation tests carried out on the parameterisation method. The focus is on the convergence of the RSVS parameterisation: first in terms of the optimality of the r-snake as a solution to the discrete length minimisation; and secondly the geometric convergence of the optimum r-snake onto the analytical solution to the governing equations.
The geometric and topological flexibility of this parameterisation is demonstrated through Fig. 7 for an aerofoil with flap. It was found that the number of snaking steps scales slowly with the number of volume cells (see Table 1 ). This is explained by the properties of the Hessian of the profile length and the Jacobian of the area constraints. Both have few off diagonal terms which means that interaction between snaxels and with the constraints is limited to those in close proximity. These properties mean that the algorithm is scalable and can be used to represent complex geometries with large numbers of snaxels and constraints. This is important as it ensures that the algorithm is capable of converging in few iterations for larger sets of design variables that could be necessary for complex topological optimisation cases.
The key consideration in allowing good performance of the algorithm is the relation between the design grid (carrying the VOS information) and the snaking grid (over which the r-snake evolves). If the underlying grid is too coarse the optimisation process cannot converge as the combination of volume constraints and smoothness conditions makes for a very stiff system. If the underlying grid is too fine there is a significant computational cost increase. For most applications a cell refinement level of 4 (each volume cell is split into 16) yields good results. Where high curvature is required within a single volume cell it can be desirable to increase the re- Fig. 7. Multi-body airfoil with parameter (volume) grid and background (snaking) grid as well as snake convergence history of the r-snake volume error and snaxel velocity. finement level, however in most cases it is preferable to increase the number of volume cells as these afford more design flexibility. However the resolution of the snaking grid controls the convergence of the RSVS on the circular arcs derived in Section 3 . Due to the discrete nature of the r-snake, the profiles generated by the RSVS only approach the shape predicted by the analytical calculations. As the number of snaxels per VOS cell is increased by using finer snaking grids the distance between the circular patches generated by NURBS and the r-snake converge. This error convergence is shown for the representation of a NACA 0012 airfoil going from a very coarse snaking grid of 2 snaxels per design variables to a finest grid of 300 snaxels per design variable in Fig. 8 . The NURBS representation is built by keeping the snaxels that lie on the edges of the design grid. These are then linked by a circular NURBS patch which satisfies the VOS requirement in the cell.
One of the key benefits of this parameterisation method over previous VOS methods is its natural ability to build sharp corners. As the required volume fraction is decreased at the edge of the profile, the minimisation of the profile length tends to create a very small feature which tends to a sharp corner as the volume fraction tends to 0. This effect can be seen in Fig. 8 , where both the leading and trailing edge are fixed in place by VOS values of 10 −5 at the extremeties. This use of small volume fractions to modify the properties of the curves is analogous tot the introduction of a knot inside a spline. A smoother leading edge can be achieved by removing these volume fractions or by designing grids with more control at the leading edge as in Fig. 9 . This ability to transition between sharp and smooth shapes is very important to the design of useful aerodynamic bodies which often require sharp trailing edges or leading edges. These small volume fractions can also be used to fix the length of a profile by effectively pinning leading and trailing edge position.

Geometric inverse design of aerofoils
To validate the geometric flexibility of the r-snake parameterisation and its suitability for aerodynamic profile generation inverse design of aerofoils was performed. The profiles were evaluated against Kulfan's Wind Tunnel Tolerance (WTT) [9] using the same process as the one used in the review of aerodynamic pa- Fig. 8. Change in normal distance between analytically derived NURBS representation and r-snake with increasing snaxel density.  rameterisation methods by Masters et al. [42] . The volume fraction values are specified on the RSVS grid by superimposing the profiles onto the grid and working out the intersection with each cell. Kulfan's WTT prescribes bounds in the maximum distance between profiles, given in Eq. (7) .
In order to achieve this tolerance for aerofoils with a reasonable number of design variables a highly anisotropic VOS grid was used. A longitudinal distribution with cells clustered at the leading edge and at the trailing edge was devised to enable accurate positioning and definition of the leading and trailing edges. In the transverse direction cells were clustered close to the chord line to allow the lower surface to cross over without causing large interferences with the upper surface. The VOS grid with the design variable values and the r-snake is shown in Fig. 9 for the case of a NACA4412 aerofoil meeting the WTT.
A detailed study of the inverse design using this VOS grid was performed on 4 aerofoils representative of common aerodynamic sections; these are the NACA 0012, the NACA 4412, the ONERA D aerofoil and the RAE 2822. The error values for these aerofoils are presented in Table 2 . Kulfan's WTT is matched for each aerofoil using the grid in Fig. 9 . This grid was tested on a further 65 NACA aerofoils of which 63 were recovered to the WTT (97%), the averaged results for this second set are in Table 3 .
While 68 design variables is more than the 20 to 30 design variables required by the established aerodynamic parametrisations studied by Masters et al. [18] for this level of fidelity; the number of active VOS design variables that need to be controlled in an optimisation is smaller than the total number of design variables in the grid used in this case. During the optimisation process only the design variables which contain the edge of the profile are of interest, this reduces the number of active design variables from 68 to 38 in the case of the NACA 4412 presented in Fig. 9 . This design variable reduction is then coupled with an overflow method which ensures smooth transition between design variables ensuring the optimisation framework only sees smooth geometric changes.
The largest errors appeared for the RAE 2822 aerofoil at the trailing edge, this is due to the very thin and curved trailing edge; this causes both the upper and lower surface to be contained in the same VOS cell which does not allow sufficient control. This case highlights the difficulty in building knowledge about a specific aerodynamic case into a very general and flexible parametrisation method. Rather than tuning the design grid to each case individually a generalised method based on local refinement of design variables is developed in Section 5 . This hierarchical approach to parameterisation offers the possibility of the RSVS tuning itself to the requirements of a given optimisation problem.

Design variable refinement
One of the drawbacks of the RSVS parameterisation method is that a regular Cartesian VOS grid contains much less implicit information about aerodynamic problems compared to traditional parameterisation methods. This means that the RSVS, while being more general than other parameterisations, also requires careful setup of the design variable layout to tackle an optimisation problem efficiently. To alleviate this, a hierarchical approach to the design variables is developed. These approaches start an optimisation problem with few design variables, this allows large but coarse changes to the design. As this process converges additional design variables are added allowing progressively finer and smaller scale changes to the design to be added. Hierarchical approaches by Anderson and Aftosmis [37] and Masters et al. [43] have accelerated and improved convergence on complex aerodynamic optimisation problems. Similar approaches have been successfully ex-  ploited in structural topology optimisation by Kim et al. to improve the performance of agent based optimisers [44] and by Bandara et al. to build a mutli-resolution framework based on sequential shape and topology optimisation using subdivision curves [45] . The RSVS lends itself to such hierarchical approaches, the refinement of design VOS design variables is intuitive and exact. A locally adaptive hierarchical process is developed in sub-Section 5.2 . Finally, the effects of local refinement are shown on inverse design and aerodynamic optimisation cases in Sections 7 and 8 respectively.

Refinement of RSVS design grids
The process of implementing local refinement into the RSVS parameterisation relies on the accurate translation of design variables from coarse to fine fidelity levels. Exact translation of a profile from a coarse to a finer layout of design variables requires only that the coarse cell boundaries is present in the finer grid. The volume fractions on the finer grid can be calculated by overlaying the profile onto the fine grid and calculating the intersection of cells and profile. This approach is analytically exact, however in practice it can lead to slightly different profiles after refinement as the snaking grid is also refined, changing the discretization of the profile.
Calculation of containments on the new design grid can be skipped by using volume fractions defined on the snaking grid calculated during the evolution of an RSVS profile. Because in most cases the snaking grid is a 4 by 4 refinement of the design grid, exact profile translation can be done using information from the snaking grid for all VOS grid refinements up to 2 by 2. The different elements of this process are presented in Fig. 10 . To fully leverage the benefits of a hierarchical approach the algorithm needs to identify and refine regions of the design space where refinement will lead to an improved objective function.

Criterion for refinement of design variables
Both uniform and local refinement of design variables have been considered by previous studies. Uniform approaches split all design variables regardless of their influence on the design; local methods aim to identify regions of the design which are more important to the reduction of the objective function and refine only those locations. Global refinement approaches have been very successful: the hierarchical approach based on subdivision curves presented by Masters et al. [43] show the best published results for the ADODG NACA0012 Case 1 [46] . Results for local refinement have been mixed, with previous studies showing that improved performance could be achieved, but the increased complexity of the optimisation method had a negative impact on the robustness of the process [43] .
Previous parameterisations that were used with refinement were tailored specifically to the aerofoil optimisation problem; this meant that the need for local adaptation by the parameterisation was limited, and that global refinement was sufficient. The goal of the method presented in this paper is to be able to tackle any shape optimisation problem with minimal tuning of the design variables by the user. For this general approach to be possible and efficient, the process developed must be able to adjust the local fidelity as the optimisation progresses. This was achieved using a local refinement algorithm which selects which VOS cells should be refined and the direction in which they should be split. The most effective criterion to be investigated refined cells with the most curvature of the profile. Curvature reliably indicates the difficulty the parameterisation is having to represent a given geometry because the minimum length objective tends to create the straightest line possible. Unlike previous methods [38,43,47] which relied on adjoint sensitivities of the objective function, this criterion relies only on information provided by the shape parameterisation.
Eq. (8) shows the calculation of the refinement criterion ( r j ) for cell j . Snaxel indices i s and i e refer to the first and last snaxel contained in the VOS cell. One critical requirement for a refinement criterion is that its value must decrease as refinement is carried out, otherwise refinement can go into a self reinforcing loop which is unlikely to lead to desirable design variable layouts. It is for this reason that coefficient r j depends on the sum of the curvature of the snake inside a cell multiplied by its length. The curvature is approximated by the second difference of the profile. The scaling using the length is needed to counteract the increasing value of curvature as feature size decreases.
with : Once a cell has been marked for refinement, the orientation of the split must be decided. Selecting the refinement direction relies on a heuristic method which yields fine grids which still enable easy movement of the profile and high quality modes. Because of the cell-bounded volume of solid formulation transition of the profile from one cell to a previously empty cell can be discontinuous. In order to minimize such transitions new cell boundaries should be normal to the profile. To achieve this, cells are split in half in the same direction as each VOS cell edge which is crossed by the r-snake. This leads to three possible refinement outcomes which are presented in Fig. 11 .

Optimisation framework
This section presents the optimisation framework that was used to test and validate the RSVS parameterisation developed in earlier sections. To test the topological flexibility of this new parameterisation each element of the optimisation framework must be able to handle the increased flexibility of the design space. The framework includes two optimisers: conjugate gradient for local optimisation and differential evolution for problems where multimodality (multiple local minima) is expected. This is coupled with an unstructured Eulerian flow solver with a cut-cell mesh generator for aerodynamic analysis.

Optimisers
The RSVS parameterisation is integrated in an aerodynamic optimisation framework which supports gradient-based and agentbased optimisers. This optimisation framework uses the volumes of the VOS cells as design variables which are translated into geometry changes by the RSVS process.
CG optimisers are easy to implement and provide adequate convergence behaviour in aerodynamic applications. The Polak-Ribiere [48] conjugate gradient formulation was adopted as it displayed better performance in early tests. The line-search is performed using backtracking until a minimum is bracketed, that minimum is then used as the step length. Differential evolution is a heuristic global optimisation method proposed by Storn and Price [49] , and was selected due to its robustness and ease of implementation both in serial and parallel. Unlike other heuristic methods it requires few internal parameters and has shown good results on a range of applications [50] , notably for constrained optimization cases.
Constraints on the design variables are handled by forcing the optimisation method to keep them satisfied at all times. For inequality constraints in addition to a hard-stop imposed on the design variables, gradients which are pointing the optimisation towards a constraint are scaled by a factor tending to 0 as the design variables tend to the constraint. This allows infeasible directions to be removed from the gradient vector and avoids the optimisation stalling close to constraints.

Smoothing of the RSVS geometric response
Desirable characteristics of parameterisation methods have been identified by previous studies into aerodynamic optimisation [18,19] . These include: smoothness, appropriate scaling and compactness. To ensure gradient based optimisers could use the RSVS efficiently the response of the parameterisation to small changes in VOS is smoothed. The smoothing removes oscillations that were found to arise from small perturbations of RSVS profiles 12 (a). This smoothing of the response is done by smearing volume fraction increments to neighbouring cells following the connectivity of the profile, this process is shown in Fig. 12 (b).
A VOS response of 25% in the first neighbours to a change in design variables was found to smooth out the oscillations on Cartesian grids. On anisotropic grids additional work is required to obtain consistent mode shapes and mode size. To smooth out oscillations on anisotropic grids the value of 25% in each neighbour is scaled by using Eq. (9) . This equation is used to reflect the different size the responses will have in neighbouring cells with very different aspect ratios and size. Finally, for anisotropic grids, the resulting smooth basis functions are scaled to have either equal height or volume depending on the objective function. This is to ensure that the resulting sensitivities are appropriately scaled.
In each VOS cell j the grid adapted coefficient ( q j ) for a smooth response is dependant on the length of the profile in the cell ( l j ) and the volume of the cell ( V j ) relative to the value for the central cell (0). It is also dependent on δ a , the size of the VOS disturbance in the central cell.

Objective functions
To test the performance of the parametrisation method on different optimisation cases the framework includes two objective functions. The first is a geometric error function which is used to perform geometric inverse design optimisations, the second is  Fig. 13. Geometric recovery of a NACA4412 using 6 refinement steps; with the RSVS grid and the target profile (left), the VOS values for the geometry (centre) and the corresponding profile coloured according to its normal distance to the target profile (right) at the first and final iteration. an inviscid Eulerian flow solver. The geometric error objective provides a fast and relatively inexpensive test case to evaluate the capability of the parameterisation method. This allows rapid evaluation of changes to the parameterisation without running a full set of aerodynamic test cases.
The geometric error is calculated as the area of the polygons created by intersecting a test geometry and a target profile. This error calculation allows the intersection of curves with different surface point distributions. This flexibility simplifies the handling of changing topology compared to an approach based on point to point distances.
Aerodynamic optimisation relies on the successful integration of CFD, parameterisation and an optimiser into a cohesive framework. For this purpose, the RSVS parameterisation method was coupled with a surface-exact cut-cell mesh generator, an unstructured Eulerian flow solver and an optimiser. In order to exploit the topological flexibility of the parameterisation all elements of the optimisation method need to support profiles made of an arbitrary number of bodies. Cut-cell mesh generators provide the required flexibility with sufficient accuracy at a low computational cost [51,52] . The flow solver is an inviscid, compressible unstructured code based on the cell-centred approach of Jameson [53] and following the implementation of Eliasson [54] . The cut-cell mesh generator and flow solver were used in previous studies by Hall et al. [55] . A mesh convergence study was performed on the zero-lift drag coefficient value for NACA 0012 at a Mach number of 0.85, giving 469.1 drag counts, which is within 0.3 counts of previous studies using different solvers [56,57] .
Because the RSVS uses traditional boundary fitted meshes it can be used with RANS solvers or any other physical solvers which uses that type of mesh. The main challenge to using the RSVS with viscous CFD is the generation of a suitable mesh without apriori knowledge of the topology. However unsupervised automatic mesh generation for viscous layers has been an active area of research recently seeing implementation in industrial codes [58] .
In the current optimisation framework the gradients are obtained by central difference on the flow solution. While adjoint gradients could also be used, the ease of implementation and parallelisation of finite differences made it a suitable option for the test cases considered in this paper. Mesh motion was carried out using the multi-scale RBF algorithm of Kedward et al. [59] as it allows efficient and exact movement of large meshes.

Geometric inverse design with refinement
The refinement process was tested on the geometric inverse design of a NACA 4412 aerofoil, using the smoothing described earlier combined with the conjugate gradient optimiser. The goal of this test is to explore the flexibility of the method and the quality of the integration with the optimiser. This case was tackled with 6 refinement steps starting from a coarse grid of 2 by 6 design variables, the grid at the first and final refinement step is shown in Fig. 13 along with the corresponding profile and volume fractions. The set-up of this case was done to test the effectiveness of the RSVS parameterisation with local refinement rather than the capacity of the RSVS to recover a NACA 4412. For this reason the profile was allowed to evolve freely over the grid with no constraint on the position of leading edge and trailing edges. The final profile ( Fig. 14 (b)) shows the parameterisation successfully built a smooth leading edge and sharp trailing edge. Building sharp trailing edges is straightforward in the RSVS parameterisation, it simply needs very small volume fractions in design cells. However this requires a design grid intersection very close to the desired location of the trailing edge which requires many refinement steps to converge as shown in Fig. 13 (b).
The same process with the same starting condition but with eight refinement stages was used to tackle the geometric recovery of a multi-body aerofoil composed of 2 NACA 4 digit profiles. Fig. 15 shows the evolution of the RSVS profile and grid through the optimisation; Fig. 16(b) the final profile with the normal dis- tance to the target profile. This second case highlights the versatility of the RSVS with refinement, their combination allows two optimisation cases with very different solutions to be tackled without user intervention.
Figs. 14 (a) and 16(a) display the convergence behaviour of the optimisations which exhibits a step by step convergence. This is the desired behaviour: as the optimisation converges on a coarse set of design variables, the refinement process selects an appropriate portion of the design variables to refine enabling further reduction of the objective function. The effectiveness of the contour curvature condition ( Eq. (8) ) is shown by the improvements brought by each refinement level; ineffective selection would lead to stagnation of the objective function.
While aerodynamic optimisation benefits from fine geometric control, it is also dependant on global parameters, especially in three dimensions global transformations such as angle of attack, sweep, twist and span must be handled concurrently to the finest design variables. In the RSVS parameterisation this can be achieved by applying those transformations to the grid so that they are reflected in the profile.
Contrary to the previous local refinement scheme developed by Masters et al. [43] which relied on the adjoint of the objective, this refinement process relies exclusively on information from the parameterisation and the profiles it generates. Instead, the objective function passes information indirectly to the refinement process through the optimisation process which drives the shape of the profile. This effect, inherent to any optimisation framework, is coupled with a measure of the 'stress' that is experienced by the parameterisation (here the curvature of the profile) to identify areas where the current parameterisation is not fine enough. This means this approach can be directly applied to any problem tackled using the RSVS parameterisation method with no modification.

NACA 0012 under local thickness contraints (ADODG Case 1)
To benchmark the optimisation framework and validate the RSVS method the ADODG case 1 was modelled [46] . The formulation of this optimisation problem is presented in Eq. (10) , where the constraint is a localised thickness constraint at every point Work in the ASO community has highlighted many of the challenges associated with the aerodynamics of this case. These difficulties include: premature convergence on sub-optimal profiles [60] , assymetric flow solutions for symmetric profiles [60,61] , oscillatory CFD solutions [62] and hysteretic behaviours with Mach number [63] . The range of observed behaviours was reviewed by Destarac et al. [64] .
This case was tackled with the CG optimiser with the mesh resolution at the 14 th refinement level which corresponds to a cell height of 0.153% of chord. This is equivalent to approximately 1300 cells uniformly distributed around the aerofoil profile. To avoid poor optimiser convergence due to assymetric flow solutions on symetric profiles, the optimisation was run on a half mesh with a symmetric boundary. Mesh-motion was performed using the multi-scale radial basis function method developed by Kedward et al. [59] . Constraints were enforced directly on the volume fractions, constraining them to be larger than required for the inverse design of the NACA 0012 aerofoil. While this does lead to a slight constraint violation in front of the point of maximum thickness it is sufficiently precise to capture the complexity of the problem. A 12 by 2 VOS grid was used: half the VOS cells are distributed in a half cosine distribution from leading edge to the 30% chord and until the trailing edge. Symmetric profiles are generated by mirroring the VOS values along the horizontal axis, meaning that the optimiser controls 12 effective design variables.
This framework allowed a drag reduction from 469 counts for the NACA 0012 to 58.3 counts for the optimised profile ( Table 4 ). This drag reduction is close to the drag values between 50 and 25 counts achieved by other aerodynamic parameterisation methods in recent comparative studies [19,64] . The shock pattern at the trailing edge of the RSVS optimised airfoil displays a supersonic/supersonic wave with a single supersonic region over the aerofoil ( Fig. 17 ). This shock pattern is similar to that observed in the review of this case performed by Destarac et al. [64] . While the drag is not as low as some previous available results, this optimisation case shows that the combination of the RSVS method with the cut-cell mesh generator is capable of exploring a complex aerodynamic design space. Flow features expected in this optimisation case are successfully discovered by the optimisation framework.

Drag minimisation for fixed area aerofoils at Mach 2
The following drag minimisations in this paper are inviscid, supersonic, constant area optimisations. Significant research into these cases was carried out in the 1950s using linearised equations for supersonic flow which yielded analytically optimal solutions. In 3 dimensions, this effort led to the now famous Sears-Haack profile for minimum wave drag [65,66] . Similar research by Klunker and Harder [67] used non-linear supersonic pressure coefficient re-lationships to obtain the profile for minimum pressure drag under thickness and volume constraints. The availability of analytical results for these cases provides useful benchmarks for non-linear numerical optimisation frameworks.
Supersonic flows are also an excellent test bed for topology optimisation: there exist multiplane profiles where shock interractions produce bodies with no wave drag [68] . The most well known of these is the Busemann biplane first proposed in the 1930s by Busemann [69] . These cases are of particular interest as these multi-body profiles can be built using the VOS parameterisation method: they are known cases for which topological flexibility should bring significant drag reduction. An example of the flow around each of the three known analytical optima is shown in Fig. 18 .
The mathematical programming representation of these problems are expressed in Eqs. (11) and (12) , for the single topology case and the multi-body case respectively. The behaviour of the optimisation is dependant on the area constraint value c A . An additional constraint is added for the multi-body cases to ensure that the optimised profile fits inside the region occupied by the Busemann biplane, the maximum height of the profile ( y max ) cannot be larger than the maximum height of a buseman biplane ( y BUSEMANN ). Indeed if the optimisation is allowed to generate profiles far apart it can effectively maintain bodies operating in flows sufficiently separated to be independent of each other.
The volume constraint was applied before the parameterisation stage by controlling the values of the volume fractions: if a constraint violations is detected the volume fractions are scaled such that their sum matches the constraint.

Impact of refinement on aerodynamic optimisation case
The first cases used to validate the current optimisation framework were the drag minimisation of profiles at Mach 2 for a fixed volume and unit chord presented in Eq. (11) using the conjugate gradient optimiser. Using the smoothed design variables the singlebody supersonic aerodynamic case was tackled for constraint values between 0.01 and 0.16 for two setups of the parameterisation.
The first was using a 2 by 10 layout of VOS design variables in a cosine distribution with symmetry of the design variables about the horizontal axis. The chord is fixed by maintaining a small volume fraction in the volume cells at the leading edge and the trailing edge. This setup is similar to traditional aerofoil parametri-sation methods where more control points are clustered towards the leading edge and trailing edges and movements in the vertical direction dominate. This approach uses engineering knowledge to build a suitable grid for the RSVS parameterisation to perform efficiently on the given problem.
The same suite of cases was then tackled using the refinement criterion specified in Section 5.2 . The optimisation was started from a 2 by 4 Cartesian grid of RSVS design variables with symmetry and five refinement steps were carried out. The goal of the refinement is to do away with the need for expert knowledge when setting up an optimisation method for a specific case. By adapting itself to the optimisation problem as it is being solved the RSVS removes a layer of complexity and improves the robustness of the entire process.
Drag results and optimum profiles for these aerodynamic cases are shown in Fig. 19 . Fig. 19 (a) shows the evolution of the drag coefficient for the analytical and non-linear optima with refinement. The inset shows the behaviour of the different optima between 0.04 and 0.048, importantly the optimisation framework successfully negotiates this complex region where two theoretical optima exist. Fig. 19 (b) shows the difference in drag value between the best analytical optima and each stage of the refinement process.
For low values of area the first refinement stages are sufficient to exceed the analytical optimum. Fig. 19 (c) shows the optimum profiles for each value of the volume constraint. The cases up to areas of 0.07 result in profiles close to parabolic but as the required volume increases, the point of maximum thickness is shifted towards the trailing edge. This finding is similar to a previous study by Palaniappan and Jameson [70] . This behaviour allows the shock to be weaker for the non-linear optimum than for the corresponding ogive, which more than makes up for the increase in back pressure. This simple geometric behaviour is easily captured by RSVS geometries with or without refinement.
For values of area above 0.09 it is very clear that the profiles tend to the truncated ogives of Klunker and Harder [67] . The large discontinuity poses a challenge to the parameterisation leading to difficulties for the optimisations. However, refinement enables the blunt trailing edge to be represented to a sufficient level and significantly improves the optimum that could be recovered compared to the cosine grid with smooth design variables. The behaviour through the refinement stages seen in Fig. 19 (b) is similar to that observed for the inverse design case presented earlier: each refinement stage unlocks a new portion of the design space to significantly improve the objective function. This helps to validate the use of curvature as a measure of the need for finer parameterisation. While there is no change in the number of bodies these results are enabled by the flexibility of the parameterisation and flow solver, traditional parameterisations are not always capable of transitioning between smooth and sharp corners as is required by the larger area cases.
The higher area cases, notably 0.15 and 0.16 exhibit a small oscillation of the profile at the trailing edge ( Fig. 19 (c)). These oscillations are the result of the optimiser minimising the turning circle of the flow to favour inviscid separation when the design variable resolution is insufficient to represent a blunt trailing edge. The large area of recirculatory flow behind the blunt trailing edge leads to poor flow convergence and poor quality gradients, preventing the optimiser from recovering at higher refinement levels. This highlights a limitation of using an Eulerian flow solver for this case: physical modelling of the boundary layer and separation is required.
The RSVS parameterisation with refinement allows the exploration of the flow for profiles of different areas in supersonic flow. The parameterisation shows that there are two different flow regimes that result in optimum solutions: attached flows for low areas and detached flows at the trailing edge for higher areas. This can be seen clearly in Fig. 19 (c): the points of maximum thickness show clear separate trends for the parabolic profiles (up to an area of 0.08) and the truncated ogives (after an area of 0.09). These results also highlight the existence of possible multi-modality around an area of 0.08.

Impact of topological flexibility on aerofoil optimisation
While the topological flexibility is not expressly needed to solve this suite of shape optimisation, it can still reveal interesting designs within the constraints of the aerofoil design space. Fig. 20 shows the flow around the best profile at iteration 4 of a fixed grid optimisation. At this point the CG optimiser has found a two body profile to be beneficial. This profile is reminiscent of the aerospike configuration used to reduce drag on the Lockheed Martin Trident D-5 submarine launched ballistic missile. This shows the capability of the topological optimisation framework as a tool for exploratory design studies. The convergence history ( Fig. 20 (a)) shows that topology change creates additional complexity in a small design space with the behaviour between iterations 2 and 6, close to the topology change, resembling convergence on a local optimum. This means that effective topology optimisation requires a method for avoiding local minima and warrants the use of a global optimiser.

Topological aerodynamic optimisation results
The approach used for the supersonic single body area constrained optimisation cases was repeated for the multi-body case. In those cases the non-linear optimum is compared to the drag value of the Busemann bi-plane. The topological flexibility is en- abled by using a larger RSVS grid layout, in this case a 10 by 6 VOS grid was used. Symmetry of the profile was enforced by mirroring the VOS design variable values: meaning the optimiser controls 30 effective design variables. To ensure exploration of the layout is effective the differential evolution (DE) optimiser is used. The starting population is of critical importance to ensure that exploration is sufficient but convergence is quick. The optimisation starts from a family of random multi-plane profiles with sharp trailing and leading edges. Starting from a population of bodies with good aerodynamic qualities reduces the convergence time significantly while still allowing the design space to be explored. A population of 100 was used as smaller populations showed inconsistant behaviour in repeated runs. The optimisation was stopped once the population showed no topological diversity. This was assessed through the convergence of the population on a set of nonzero VOS design variables, This occurred between iterations 200 and 300. Fig. 21 presents the drag results and profiles resulting from the topological optimisation process. Fig. 21 (b) shows the drag of the optimised profiles is below that of the Busemann biplane or the best analytical optima for all values of the constraint above 0.02. This good performance above 0.02 is because the optimisation tends to build very efficient profiles which resemble convergent divergent nozzles with flat outer edges. The smooth compression which results from these has a much lower drag on a discrete grid compared to the Busemann bi-plane which relies on perfect shocks and expansion fans. Below 0.02 the profiles are extremely thin, building planes less than 1% thick. This means that the sharpness of the leading and trailing edge play an outsize role in the quality of the optimum, however the optimiser struggles to adjust these sufficiently.
Above an area of 0.1 the flow in the Busemann bi-plane is choked which can be seen in the very large increase in drag (18 times larger between areas of 0.1 and 0.11). In Fig. 21 (c) the profile of area 0.12 displays the internal features similar to the optima of area 0.0315 but with curved outer edges. Fig. 22 (c) shows that the optimiser is combining the flow features of an optimised Busemann bi-plane (flow similar to Fig. 22 (a)) with one of the single body optima of Fig. 19 (flow similar to 18 (c)). This allows low drag to be maintained where a traditional Busemann bi-plane would choke.
The optimum profile for an area of 0.06 shows a pentaplane profile ( Fig. 22 (b)), while for 0.0315 the optimum is a tri-plane ( Fig. 22 (a)). This difference in optimum topology is because the main factor in drag reduction is the minimisation of the external shocks which can be achieved with either topology. This also explains the large number of iterations required to converge the topology of the optimisation case in Fig. 23 : the optimum topology is not stable before the 234 th iteration.
These differences in optima highlight a limitation of the differential evolution on this topological aerodynamic optimisation case. While DE provides good exploration, the convergence on the global optimum is not guaranteed; this is because each of the local minima has a very similar drag value but with different topology. Alternate algorithms for niching and hybrid gradient/agent search methods could help improve the performance of the framework on these cases. Despite these limitations, the combination of parameterisation, global optimiser and flow solver is effective at exploring these optimisation problems. The relative compactness of the set of design variables as well as the smoothness of the recovered profiles ensures that good aerodynamic bodies are generated most of the time without arbitrarily restricting the design space. This allows this topological optimisation framework to explore the complex behaviour of the optimal solution for large values of area.

Conclusions
Optimisation for external aerodynamics has usually focussed on small surface changes, but including topological change within the calculation makes accessible a rich landscape alternative designs. In some circumstances, exploring these widely varying alternatives thoroughly can be critical for the success of the optimisation.
To achieve this flexibility, a technique using a local fraction of solid across a design parameter grid has been developed for representing two-dimensional shapes, including changes in topology. Although a coarse parameter grid may be used initially, sequential anisotropic refinement of the parameter grid is able to reconstruct aerofoils to a good degree of accuracy. The method does not achieve the same level of efficiency (in terms of numbers of design variables required) compared to parameterisation methods specific to aerofoils, but it does create a flexible, general method for any shape.
The construction of the surface is implemented using r-snakes, and the snaxel positions of the r-snakes are determined so as to match each specified cell volume of solid while also minimising the surface arc length. Calculus of variations shows that this for-   Fig. 21 (c). mulation leads to a surface built of circular arcs, and therefore the surface may also be represented exactly using NURBS if required. Small perturbations in the parameters lead to parabolic shape changes, and these permit a system of smoothing to be constructed.
When tested for geometric reconstruction, the method is able to reconstruct both single and multiple aerofoil sections. Fixed volume drag minimisation in supersonic flow has revealed known parabolic and truncated shapes for single shapes, and when multiple shapes are allowed, Busemann multi-planes are generated. Fu- ture work will extend the technique to three dimensions and explore the potential for other search methods to explore aerodynamic topology optimisation cases more effectively.

Acknowledgements
This work was carried out using the computational facilities of the Advanced Computing Research Centre, University of Bristol -http://www.bris.ac.uk/acrc/ . Funding: This work was supported by the EPSRC Research Council UK and the University of Bristol through DTP grant EP/M507994/1 . All the data needed to support the conclusions of this article are presented in the paper.

Appendix A. Calculation and differentiation of the volume fraction
Calculation of the volume fraction is performed using Green's theorem which allows the calculation of the area of closed curves. The two dimensional form of the theorem for a polygon is presented in Eq. (13) ; the results of this theorem can be extended to higher dimensions for a 3 dimensional implementation of this method. q is the centre of the edge of polygons and n the outward facing normal to the edge.
This equation is used within each cell to define the area currently contained by the snakes. Manipulation of Eq. (13) allows it to be formulated into the matrix product in Eq. (14) . Vector p is the list of coordinates of the vertices where p n is the coordinate row vector of the n th vertex defining polygon C . Eq. (14) is derived from the decomposition of q and n into, respectively, the mean and the difference of neighbouring vertices which is readily transformed into simple matrix equations that can be assembled to reflect the connectivity information of a polygon. More detail of the derivation of matrix R A is presented in Appendix B .
The SQP algorithm being developed requires the derivation of the Jacobian of the area constraints with respect to the nondimensional snaxel distances. Applying the matrix form of Green's theorem to the VOS cells greatly simplifies the calculation of the derivatives as all the connectivity information is precomputed and hidden into matrix R A . For the derivative to be computed, vector p is readily separated into a variable and constant part. Recalling the formulation of p i ( Eq. (2) ) and of p ( Eq. (14) ). 1 , g 2 , 1 , · · · , g n, 1 } T (15) This formulation shows that p is only a function of the distance along edges of the snaxels ( d i ) and properties of the snaxel grid, the direction of a travelling snaxel ( g i ) and the originating grid point ( g i ,1 ). For any vertex of the polygon which is not a snaxel the entry into d C is replaced by a 0. Replacing Eq. (15) into Eqs. (14) and (16) is developed. This form of the equation simplifies the differentiation process with regard to d .
The differentiation relies on simple matrix derivation rules and the symmetric nature of [ R A ]. In the previous equations I d is the result of the operation ∇ d d C ; it is a rectangular matrix of ones and zeros which has the effect of deleting rows and columns from the equations corresponding to static vertices and inactive snaxels.
∇ d A C is the gradient of a VOS cell; it is a column vector of length n (the number of snaxels in the profile). To build the full Jacobian of the constraint ( ∇ d h ) the gradient in each VOS cell is calculated using Eq. (17) and the resulting vectors are assembled to form the matrix of Eq. (19) .  with : This form allows a much more readable representation of the first and second derivatives of the function.
is a small positive number used to stabilise the derivatives as the edge length goes to zero. Only the properties of the Hessian and Jacobian are discussed in this section; the equations for the derivatives of Eq. (23) are presented in Appendix D . The Hessian is a tridiagonal symmetric matrix ( Eq. (18) ) which means the cost of inverting it for the calculation of the Newton step ( Eq. (3) ) is low. This is due to the formulation of the tensile force, as it only relies on one neighbour on each side it leads to a sparse Hessian which favours the stability of the system. The value of ( Eq. (23) ) is chosen to ensure that the denominator of the derivatives does not go to 0 and is sufficiently high to ensure good conditioning of the Hessian. A value is selected such that the impact on the derivative is limited to a small region around singularities. A typical value for this parameter will be of the order 10 −5 .

Supplementary material
Supplementary material associated with this article can be found, in the online version, at doi: 10.1016/j.compfluid.2019.02. 008 .