INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN FLUIDS Int. J. Numer. Meth. Fluids (2013) Published online in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/fld.3844

CAD-based shape optimisation with CFD using a discrete adjoint Shenren Xu1 , Wolfram Jahn2 and Jens-Dominik Müller1, * ,† 1 School

of Engineering and Materials Science, Queen Mary, University of London, London, UK 2 BDSP Partnership, London, UK

SUMMARY One of the major challenges of shape optimisation in practical industrial cases is to generically parametrise the wide range of complex shapes. A novel approach is presented, which takes CAD descriptions as input and produces the optimal shape in CAD form using the control points of the Non-Uniform Rational B-Splines (NURBS) boundary representation as design variables. An implementation of the NURBS equations in source allows to include the CAD-based shape deformation inside the design loop and evaluate its sensitivities efficiently and robustly. In order to maintain or establish the required level of geometric continuity across patch interfaces, geometric constraints are imposed on the control point displacements. The paper discusses the discrete adjoint flow solver used and the computation of the complete sensitivities of the design loop by differentiating all components using automatic differentiation tools. The resulting rich but smooth deformation space is demonstrated on the optimisation of a vehicle climate duct. Copyright © 2013 John Wiley & Sons, Ltd. Received 18 February 2013; Revised 20 July 2013; Accepted 30 August 2013

KEY WORDS:

shape optimisation; discrete adjoint; CAD; NURBS

1. INTRODUCTION Numerical shape optimisation has been applied in aeronautical engineering for many decades [1, 2] and is now maturing in a wider range of applications, such as the automotive industry [3, 4]. The currently used approaches fall into two main categories: stochastic and gradient-based methods. Stochastic methods, such as genetic algorithm and evolutionary methods [5], are robust and easy to implement using black-box components but have prohibitive computational cost when applied to practical industrial cases with many design variables. The gradient-based approach, on the other hand, can be more efficient, especially if there are many design variables and if the baseline design is already close to optimal, both of which are typically the case in industrial applications. The main challenge in gradient-based methods is to evaluate the gradient in a computationally efficient manner. The adjoint approach [1] has become the method of choice for computing gradients as the computational cost of gradient computation becomes virtually independent of the number of design variables. The significant disadvantage of the adjoint approach is that in addition to the CFD solver, an adjoint solver needs to be developed. Early work on adjoint development [2] proposed to derive the adjoint equations and then discretise them, the ‘differentiate then discretise’ or ‘continuous’ approach. This approach is also used in the open-source solver OpenFOAM [4]. Because of the re-discretisation of the problem, the obtained gradient is the inexact gradient of the exact cost function; that is, the location of the minimum of the cost function predicted by the flow solver will only coincide with zero gradients of the adjoint at infinitely fine mesh resolution. *Correspondence to: Jens-Dominik Müller, School of Engineering and Materials Science, Queen Mary, University of London, Mile End Road, London, UK. † E-mail: [email protected] Copyright © 2013 John Wiley & Sons, Ltd.

S. XU, W. JAHN AND J.-D. MÜLLER

More recently, work has focused on basing the adjoint solver on the exact Jacobian of the flow solver by differentiating the statements in the CFD code implementation, the ‘discretise then differentiate’ or ‘discrete’ approach. The differentiation of the flow solver, the ‘primal’, can be performed by hand [6], but advances with automatic differentiation (AD) software tools [7] can be exploited to automate this process, as the application of the rules of calculus to each statement can be applied methodically. The discrete approach computes the exact sensitivity of the inexact cost function; that is, the computed gradients always match the solution given by the flow solver. Because the problem is not re-discretised, the truncation error of the cost function in primal and adjoint is the same. Arguably, the discrete adjoint is actually more appropriately seen as the differentiation of an algorithm, which in this case approximates the solution to the flow equations, rather than a solver for the adjoint equations. Application of AD is not as simple and straightforward as suggested earlier for large codes, especially if one requires code with adequate performance in runtime and memory use. Extensive adaptation of the primal is often required [7–12]. Another major concern in industrial CFD shape optimisation is the choice of parametrisation [13, 14]. For routine industrial application, we need an approach that (a) is fully automatic, that is, does not require additional user input; (b) supports a rich design space of smooth deformations, that is, rich enough to capture all relevant modes while not permitting oscillatory modes that may be difficult to damp; (c) permits the imposition of geometric constraints, for example, due to buildspace or manufacturing; (d) returns the geometry in a CAD format back into the virtual prototyping chain for further analysis and manufacturing; and (e) can be integrated into the design loop robustly with low computational cost. A number of approaches have been proposed: node-based [15, 16], radial-basis function [17], shape functions on lattices [18], free-form deformation [19] and CAD-based parametrisation [20]. Unfortunately, none of the currently proposed methods satisfy all the aforementioned requirements. Lattice-based shape functions [18] and free-form deformations [19] require auxiliary shape grids; their construction becomes extremely cumbersome when rich design spaces are needed for complex parts. In these methods, the numerical conditioning of the matrix representing the derivatives of surface node displacements with respect to design variables is often poor as shape modes are not sufficiently independent of each other. The human effort in setting up the design space typically results in a restricted design space that may not contain relevant minima [19, 21]. In addition, the optimal design again exists as a deformed mesh, and a manual transformation back to CAD is required, which may lose some of the relevant shape modes and hence affect the optimum. Radial-basis function approaches [17] simplify setting up a smooth deformation field as no lattice is required, but it is difficult to constrain shape deformations between the basis function points. Node-based deformation uses the displacement of the nodes of the surface grid as design variables. The method is automatically set up through the surface grid and represents the richest design space that can be expressed on the CFD mesh. However, high-frequency oscillatory modes resulting from reduced regularity of the gradients [15] are in general not sufficiently resolved by the CFD solver to be visible in the cost function and are hence not appropriately suppressed by the optimiser. Additional smoothing is required [15, 16, 22]. The major drawback is that the optimised shape exists only as a deformed surface mesh; again a manual transformation back to CAD is required with the same disadvantage as for the other aforementioned methods. However, the rich design space and easy set-up makes the node-based approach a very valuable tool for design space exploration. For a complete optimisation workflow to be used in routine design, the optimised shape needs to be produced in a CAD format. A first CAD-based approach is to directly employ the parametrisation defined in a CAD system [23, 24]. These parametrisations are typically built to generate families of parts and in general will need to be substantially modified to produce a suitable design space to capture the optimum, hence still requiring significant user input. The manual set-up restricts the design space, and an important mode may not be represented. On the other hand, the provided parametrisations often have the relevant geometric constraints built in. There are however severe shortcomings with this approach. Firstly, there is no universal standard for these CAD parametrisations; they are defined using standards specific to a CAD system that are in general not transferable. Secondly, and more importantly, current CAD systems do not offer derivatives of surface displacements with Copyright © 2013 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Fluids (2013) DOI: 10.1002/fld

CAD-BASED SHAPE OPTIMISATION WITH CFD USING A DISCRETE ADJOINT

respect to the design parameters, which are needed in the chain rule to compute the sensitivities of the cost function with respect to design variables. The only available option is to apply finite differences to the CAD system, resulting in significant robustness issues as finite-size displacements often lead to topological changes in the surface description and changes in patch numbering. Robinson et al. [24] also report significant runtimes for a 3D air duct case with 174 design variables. If one uses an efficient simultaneous time-stepping method [25, 26], this penalty would increase dramatically as these methods use a large number of small design steps. We propose an alternative approach that avoids the robustness and runtime issues of using the CAD parametrisation. The approach, first presented in 2D by Yu et al. [27], utilises the fact that the relevant output of the CAD system to the analysis and manufacturing tools in the virtual prototyping chain is the boundary representation (BRep) given in the STEP standard as a set of Non-Uniform Rational B-Splines (NURBS) surface patches [28]. The approach uses the displacements of the control points of the NURBS patches as design variables; the updated set of control point positions fully define the CAD shape of the optimal design at convergence. This approach produces a modified CAD description of the surface at output. It is fully automatic as it does not require any additional user input and it works with a generic, vendor-independent parametrisation. A number of authors have used the NURBS control points as design variables [21, 29] but considered only geometries with a single patch or only deformation inside a patch. The main challenge to generalise this approach is to maintain the required level of continuity of tangency and curvature between adjacent NURBS patches when a control point on or near a patch interface is displaced. The proposed method introduces constraints for geometric continuity across NURBS patch interfaces and is hence termed ‘NURBS-based parametrisation with continuity constraints’ or NsPCC. This paper first explains the discrete adjoint solver in Section 2, followed by a description of the mesh smoothing in Section 3. Section 4 discusses the formulation of the CAD-based parametrisation. Its implementation and how to satisfy and maintain constraints of geometric continuity is discussed in Sections 5 and 6. Finally, Section 7 analyses the results of the optimisation of a Volkswagen air duct and demonstrates the effectiveness of the methodology. Conclusions are presented in Section 8. 2. THE ADJOINT SOLVER The fitness of a design can be measured with a scalar cost function J.U.˛/, ˛/, which depends on the state U , which in turn depends on some design variables ˛. In gradient-based optimisation, one needs to compute the gradient dJ =d˛ of this cost function with respect to the design variables at each design step,   n  dJ n n ˛ CA ˛ , ! ˛ nC1 , (1) d˛ where the operator A.˛, rJ / represents an optimisation algorithm that returns a perturbation ı˛ to the current design ˛ n . This process is iterated until ˛ nC1 meets a stopping criterion, such as the magnitude of the gradient being acceptably small. The minimum that is sought is subject to the steady-state PDE constraint of the flow equations, which can be written as R.U.˛/, ˛/ D 0.

(2)

Differentiating (2) with respect to a design variable ˛ results in @R @R @U C D 0, @˛ @U @˛

(3)

Lu D f ,

(4)

which can be written in short form as

where LD

@R , @U

Copyright © 2013 John Wiley & Sons, Ltd.

uD

@U @˛

and f D 

@R . @˛ Int. J. Numer. Meth. Fluids (2013) DOI: 10.1002/fld

S. XU, W. JAHN AND J.-D. MÜLLER

The gradient thus can be computed as @J dJ D C g T u, d˛ @˛

(5)

where gT D

@J . @U

The perturbation field u is obtained from solving the linear equation (4) for each design variable, as the right-hand side f depends on ˛, dJ @J D C g T L1 f . d˛ @˛

(6)

This is termed the tangent-linear approach. To derive the adjoint approach, we regroup the transpose in the field term of (6) as dJ @J  T T T D C f L g . d˛ @˛

(7)

As in the tangent-linear approach, the term LT g corresponds to solving a linear system for the adjoint variable v, LT v D g,

(8)

dJ @J D C vT f . d˛ @˛

(9)

which simplifies the adjoint sensitivity to

The equivalence between the terms g T u and v T f in the sensitivity equations (5) and (9) can be easily validated with the following identity that holds if the Jacobian L is transposed exactly:  T v T f  v T Lu  LT v u  g T u. The adjoint equations (8) depend on the cost function through the right-hand side g but are independent of the design variables ˛. The adjoint sensitivity (9) does, of course, depend on the design variable through the source term f , but f does not involve a linear system solve and v Tf can be evaluated efficiently using reverse differentiation. With the adjoint approach, the cost of computing sensitivities becomes independent of the number of design variables, which is the enabling feature for automatic parametrisation approaches that typically result in very large design spaces. In this work, we use an in-house incompressible discrete adjoint solver, gpde [12], developed using the AD tool Tapenade [9]. The primal flow solver of gpde uses a standard incompressible method with cell-centred storage, face-based residual assembly and, as an iterative scheme, the SIMPLE pressure-correction discretisation [30]. The stabilised bi-conjugate gradients (Bi-CGSTAB) method [31] is used to solve the linear system in the momentum step of the SIMPLE scheme, and the conjugate gradient (CG) method is used to solve the pressure correction equation. The adjoint solver is assembled by treating the (pseudo) time-steps of the incompressible flow solver as a fixedpoint iteration, that is, the sensitivity computation does not depend on the convergence history. As a consequence, the primal (flow) solution does not need to be recorded (checkpointed). Checkpointing is required for the outer iterations within one cycle of momentum prediction and pressure correction as the SIMPLE scheme solves for each block separately. The backward differentiation is performed by applying the AD tool to this outer iteration. The inner iterations for each block, that is, the linear solves for momentum and pressure correction, are straightforward linear solves with a perturbed right-hand side; hence, the adjoint equivalent is the same linear solve with a transposed system matrix and different right-hand side. Differentiation of the linear solver is possible Copyright © 2013 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Fluids (2013) DOI: 10.1002/fld

CAD-BASED SHAPE OPTIMISATION WITH CFD USING A DISCRETE ADJOINT

but is inefficient as the differentiation would store the basis vectors of the Krylov method used by the primal solve and apply those to solve the tangent or adjoint linear systems. A better approach is to call the same linear solver with altered arguments, basing the nonlinear steps in the Krylov solves of the adjoint system on iterates of the actual adjoint solution. Tapenade allows the user to specify a provided derivative code that is to be used instead of the differentiated code. We chose here to replace Tapenade’s differentiated linear solver code with a hand-coded one using scripting in a post-processing step. The resulting adjoint code then has a runtime between four and five times of that of the primal solver, but this performance may be further improved by code tuning [11]. 3. MESH DEFORMATION ALGORITHM At each design step, the modified deformation of the surface needs to be propagated smoothly into the domain to maintain the mesh quality at a sufficient level. Maintaining mesh quality under large deformations is a critical aspect for application to industrial cases, partly to maintain numerical accuracy, but primarily to maintain the stability of the solver throughout the optimisation. The mesh deformation used in this work is based on a modified version of the spring analogy and proves to be sufficient for the relatively large deformation in the presented optimisation case. The generic form of the spring analogy considers each edge of the mesh to be a spring with nodes in equilibrium if the sum of spring forces is zero. Taking the spring force proportional to the design displacement of the node leads to the following equation: ıXi D

kij ıXj †ngh j †ngh ki ,j j

,

(10)

where the sum runs over all ngh neighbouring nodes connected to i through an edge, kij denotes the spring stiffness between nodes i and j , and ıX is the displacement of each node with the displacement of surface nodes ıXs imposed by the optimiser. The governing equation can be written as ŒKıX D ıXs ,

(11)

which is solved efficiently using the CG method at a computational cost that is negligible compared with that of the primal and adjoint solvers. To better preserve the quality of highly anisotropic boundary layer meshes, the spring stiffness often is made proportional to the inverse of the distance between two nodes as kij D

1 , jXi  Xj j

(12)

which creates a ‘repulsive’ force by making very short edges very stiff [32]. This is further modified here to kij D

dW Amax ,  jXi  Xj j Ami n

(13)

where dW means the distance from the node to the wall as required by, for example, the Spalart– Allmaras turbulence model, and Amax and Ami n mean the largest and smallest areas among all the faces of one element. This scaling maintains the typically used inverse relation with edge length and also makes the stiffness of edges proportional to wall distance and aspect ratio evaluated using face areas. With this modification, the perturbation of the boundary will move the boundary layer mesh almost as a rigid body and propagate the displacement into the domain far from the wall. @J In discrete adjoint methods, the sensitivity @X of the cost function J with respect to a s perturbation in surface node coordinates Xs is written as @J D v Tf , @Xs Copyright © 2013 John Wiley & Sons, Ltd.

(14) Int. J. Numer. Meth. Fluids (2013) DOI: 10.1002/fld

S. XU, W. JAHN AND J.-D. MÜLLER

with vD

@R T @J @U @U

and f D

@R @X . @X @Xs

(15)

This product can be evaluated efficiently in a consistent way by first solving the adjoint equations ATv D g and then computing the product v Tf by reverse differentiating the metrics calculation @R (face/edge normals) to obtain @X and then chaining that with reverse differentiation of the same @X mesh smoothing algorithm of (13) to right multiply with @X [33]. s The metrics and mesh deformation subroutines are differentiated using the AD tool Tapenade in reverse mode for use in the adjoint code. The sensitivity of the cost function with respect to the volume grid is used as input for the reverse-differentiated mesh deformation subroutine, and the sensitivity of the cost function with respect to the design surface mesh will be the output denoted as dJ =dXs . 4. CAD-BASED PARAMETRISATION The CAD-based parametrisations evaluate the deformation of the surface resulting from a perturbation in design variables. The surface deformation then needs to be interpolated onto the surface nodes of the CFD mesh and propagated smoothly into the domain to maintain volume mesh validity and quality. The smoothing approach taken here has been discussed in Section 3. Starting from Equations (14) and (15), the adjoint-based sensitivity (9) can then be written as dJ dJ dXs D , d˛ dXs d˛ where Xs is the surface node coordinate and ˛ is a CAD parameter. s Robinson et al. [24] use finite differences to compute the shape derivative dX , where ˛ are d˛ user-defined and system-defined CAD parameters. The main issues with this approach have been discussed in Section 1. In contrast to this, we propose to base the CAD parametrisation on the BRep, which uses a collection of NURBS patches as given in the standardised STEP format. To modify the shape, each control point of a NURBS patch is allowed to move in all directions, hence each representing 3 degrees of freedom (DoF). A simple example of a NURBS patch with its associated controls points is shown in Figure 1, which also illustrates how the perturbation of a control point changes the shape of the NURBS patch. Allowing every control point to move in all directions represents the richest design space the CAD model is able to express. Adaptive control point enrichment or order elevation of NURBS patches can be implemented to refine the design in areas of high sensitivity.

Original NURBS

Perturbed NURBS

Figure 1. A NURBS patch with the net of original (upper left) and perturbed (lower right) control points. Copyright © 2013 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Fluids (2013) DOI: 10.1002/fld

CAD-BASED SHAPE OPTIMISATION WITH CFD USING A DISCRETE ADJOINT

Imposing geometric buildspace constraints can be carried out either at the level of the control points exploiting the convex-hull property of splines [28] or alternatively at the level of the surface mesh points, which achieves tighter bounds at higher computational cost. In addition to requiring smooth shapes, industrial applications often seek to limit the local curvature, for example, at trailing edges of aero-engine turbine blades. Unfortunately, there is no known method to derive analytic bounds for the maximal curvature in a NURBS patch, but control of maximal curvature can be performed at the surface mesh level, which is the appropriate level of detail resolved by the CFD simulation. These build constraints can be implemented in the same framework as the patch continuity constraints discussed in Section 5. For an independently moveable control point with coordinates P , the gradient of the cost function can then be written as dJ dJ @Xs D . (16) dP dXs @P s The shape derivatives @X can be calculated analytically; their definition is independent of the CAD @P system. Specifically, a NURBS is defined [28] as

Xs .u, v/ D

n X m X

Bi ,j .u, v/Pi ,j ,

(17)

i D0 j D0

where Pi ,j are the control points, and u and v are the parametric variables of the surface mesh point. The rational basis functions Bi ,j .u, v/ are defined as Ni ,p .u/Nj ,q .v/wi ,j Pm , kD0 lD0 Nk,p .u/Nl,q .v/wk,l

Bi ,j .u, v/ D Pn

(18)

where Ni ,p .u/ and Nj ,q .v/ are the p-th and q-th degree basis functions defined on the parametric space u, v. Note that the shape derivative matrix @Xs =@P is equal to the basis functions Bi ,j , and thus, the gradient of the cost function with respect to the control points is simply dJ dJ B, D dP dXs

(19)

where B is a matrix with elements Bi ,j . In addition to the shape derivative, we also need to compute actual displacements following a shape update by the optimiser. We have hence chosen to implement a basic CAD kernel that evaluates NURBS surfaces; this code has then been differentiated using the AD tool Tapenade [9] in forward mode to provide the necessary derivatives. 5. GEOMETRIC CONTINUITY The finite displacement of a control point P on or near a patch interface typically results in violation of the continuity constraints; for example, control points on an interface to a fixed/nonmoveable patch must not move at all to maintain G0 continuity (no gaps). Requiring in this case G1 (tangent) and G2 (curvature) continuity will additionally lock the second and third rows of control points. Similarly, control point displacements on moveable patch interfaces imply constraints on the neighbouring rows of control points along the interfaces. In general, G0 continuity requires that the adjacent patches share a common edge. G1 continuity requires firstly that G0 is satisfied and, secondly, that the adjacent NURBS patches share the same tangent plane for any point along the common edge. However, this does not necessarily require that the adjacent patches have the same number and distribution of control points or degrees of basis function. In our approach, the constraints are evaluated at a number of test points that are distributed along a coincident parametric edge or the intersection line of both patches. Figure 2 shows two NURBS patches sharing one common edge. Note that the number of control points along the common edge could be different for the left and right patches, but the test points are always deployed in pairs, with one on each NURBS patch. The required number of test point pairs along an edge is disCopyright © 2013 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Fluids (2013) DOI: 10.1002/fld

S. XU, W. JAHN AND J.-D. MÜLLER

test point

control point

Figure 2. Test points along the NURBS boundaries and the relevant control points on each patch.

cussed in Section 5.1. Continuity constraints for each pair of test points then express that location, tangent plane (for G1 ) and curvature (for G2 ) are identical on all patches containing this pair of test points. For example, for G1, the position vectors and the normal vectors of the tangent planes need to be identical at the test point when evaluated in either patch. A kink between patches, for example, as a manufacturing constraint for mould deforming, can be specified by requiring normals to differ by the imposed deforming angle. For G0 and G1 continuity along an edge shared by two NURBS patches, the following constraint functions need to remain zero for the test points evaluated in each of the patches: G0 D Xs,L  Xs,R

(20)

G1 D nE L  nE R ,

(21)

where the discretely evaluated positional and tangent continuity constraint functions are denoted as G0 and G1 , respectively; note the difference in notation between the discrete evaluation G0 of Equation (20) and the analytic form denoted as G0. The unit normal nE of the tangent plane at the test point is defined as   @Xs @Xs  @u @v ˇˇ , nE D ˇˇ (22) ˇˇ @Xs @Xs ˇˇ ˇˇ @u  @v ˇˇ where u and v are the surface parametric variables as introduced in Equation (17) and the subscripts L or R mean that the term is evaluated for the test point on either the left or right patch, with the notation suitably extended for corner vertices. To maintain the geometric continuity of the initial geometry, each constraint function is required to remain zero after each design update. To evaluate the change in constraints, each constraint function is linearised as G nC1  G n C

N X @G ıPi , @Pi

(23)

i D1

where Pi denotes the displacement of the x, y and ´ components of the i-th of N control points. Equation (23) uses the superscripts n, nC1 to express the numerical evaluation of the constraint at the design iterations n, nC1. Assembling (23) for each test point displacement, and requiring G n  G nC1 D 0, the following linear equation system is obtained: CıP D 0, Copyright © 2013 John Wiley & Sons, Ltd.

(24) Int. J. Numer. Meth. Fluids (2013) DOI: 10.1002/fld

CAD-BASED SHAPE OPTIMISATION WITH CFD USING A DISCRETE ADJOINT

where

2

@G01 @P1

6 6 6 1 6 @G1 6 @P1 6 6 6 " i # 6 @G02 6 @P1 @Gj 6 CD D6 6 @G 2 @Pk 6 1 6 @P1 6 6 6 6  6 6 4 M @G1 @P1

@G01 @P2



@G11 @P2



@G02 @P2



@G12 @P2

 ..

 @G1M @P2

.



@G01 @PN

3

7 7 7 @G11 7 @PN 7 7 7 7 2 @G0 7 @PN 7 7 7 2 7 @G1 7 @PN 7 7 7 7  7 7 7 5 M

2

and

ıP1

6 6 6 ıP2 6 ıP D 6 6 . 6 . 6 . 4

3 7 7 7 7 7. 7 7 7 5

ıPN

@G1 @PN

if G1 is imposed. In this notation for G, the superscript i refers to the i-th of M pairs of test points, and the subscript j denotes the Gj constraint function (j D 0 for positional constraint and j D 1 for tangent constraint). The subscript k to Pk refers to the displacement of the k-th of N NURBS control points. The linearised constraint matrix C has 23M rows corresponding to a total of M pairs of test points with G1 constraint and has 3N columns corresponding to a total of N control points. To satisfy the continuity constraints in a linearised sense, the perturbations of the control points ıP have to lie within the null space of the constraint matrix C. We compute the null space of C with an SVD: C D U†VT ,

(25)

where U is a mm unitary matrix, † is a mn diagonal matrix with non-negative real numbers on the diagonal and the nn unitary matrix VT denotes the transpose of V. The rank of the matrix, r, is the number of the non-zero diagonal entries in †. The last .m  r/ columns of the matrix V span the null space of C, denoted by Ker.C/. As rounding error may lead to small but non-zero singular values, it is necessary to find a demarcation point below which the singular values are regarded as zero. The SVD orders the diagonal elements of † by the magnitude of the singular value i , ranking highest those modes that are most important to approximate C. This property enables us to truncate the null space at an appropriate level as discussed in Section 5.1. The allowable control point perturbations hence become ıP D

mr X

vkCr ı˛k D Ker.C/ı˛,

(26)

kD1

where ı˛k , k D 1, : : : , m  r are the perturbations in design parameters. The gradient formulated in Equation (16) can then be further modified to be dJ dJ @Xs @P dJ @Xs D D Ker.C/. d˛ dXs @P @˛ dXs @P

(27)

5.1. Required number of test points As NURBS is based on polynomials over specific knot-vector intervals, testing the values of two NURBS at a sufficient number of points will allow to test exactly whether they match. It is however not straightforward to pre-determine the number of test points that are needed [28]. We propose to use the number of non-zero singular values of the linearised constraint matrix C to determine the number of test points needed along each edge. The SVD will filter out redundant constraints arising from excessive test points: once there are sufficient test points to determine the null space, the number of singular values will no longer increase when a further control point is added. Copyright © 2013 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Fluids (2013) DOI: 10.1002/fld

S. XU, W. JAHN AND J.-D. MÜLLER

Figure 3. Left: Singular values for the linearised G1 constraint matrix using different number of test points along the edges; right: the number of allowable displacement modes according to eigenvalues-value analysis with cut-off value of  D 107 .

Singular values are considered zero if they fall below a certain threshold. G0 constraints lead to a very cleanly conditioned SVD with singular values either clearly non-zero or of the order of machine precision. In the case of G1 or higher, the distribution of singular values is more uniform with a long tail of small singular values gradually tailing off in size. We select a cut-off for the singular value at 107 , which may reduce the accuracy of SVD but increases the design space. As G1 and higher constraints are nonlinear, the design step in the null space will not satisfy the constraints exactly, and a recovery step is needed, as shown in Section 6. A very exact representation of the null space is hence not required. This methodology is illustrated using the S-Bend optimisation case described in Section 7. When only G0 is imposed, the number of non-zero eigenvalues remains unchanged if there are more than 12 test points (figure not shown here). Figure 3 shows the behaviour for G1; the number of non-zero singular values ( > 107 ) ceases to increase for more than 30 pairs of test points along each edge. 6. G1 CONTINUITY RECOVERY The geometric continuity constraint developed in Section 5 restricts the perturbation to be tangent to the linearised constraint functions. The G1 constraint function is nonlinear; hence, each linearised tangent step will slightly violate the constraint. To recover the constraint after each tangent step, additional normal steps in the range space of C are required to recover the nonlinear continuity constraints G1 and G2. The following development shows recovery of G1, and the extension to G2 is straightforward. Defining the deviation of the G1 constraint function from the target by ıG1 ,

ıG1 D

M X

jjG1i jj,

(28)

i D1

where M is the number of pairs of test points and Gi1 is the constraint function evaluated at the i-th pair of test points, and denoting the linearised G1 constraint function by C1 as in Equation (24), we can compute the recovery step of the control points ıP? as C1 ıP? C ıG1 D 0. Copyright © 2013 John Wiley & Sons, Ltd.

(29) Int. J. Numer. Meth. Fluids (2013) DOI: 10.1002/fld

CAD-BASED SHAPE OPTIMISATION WITH CFD USING A DISCRETE ADJOINT 0.06 Without G1 recovery With G1 recovery

Deviation from G1

0.05 0.04 0.03 0.02 0.01 0 -0.01

0

2

4

6

8

10

Optimisation iterations

Figure 4. Deviation from exact G1 continuity with versus without G1 recovery steps.

In addition to minimising G1 , we also require G0 to be still strictly satisfied, and thus, ıP should satisfy during each recovery step ıP? D Ker.C0 /ı˛? ,

(30)

0 is calculated and stored at the beginning of the whole optimisation, as C0 is where C0 D @G @P independent of control point displacements. Algorithmically, we first solve for ı˛? , which satisfies

C1 Ker.C0 /ı˛? C ıG1 D 0,

(31)

using SVD for pseudo inversion of the system matrix C1 Ker.C0 /. As a second step, we compute the control point corrections ıP? using Equation (30), which guarantees that they lie strictly in the null space of C0 , and thus, G0 continuity is exactly satisfied in the G1 recovery step. In the test case presented in this work, only a few Newton steps are needed to bring the G1 continuity below an appropriately chosen threshold. The effect of applying G1 continuity recovery steps is shown in Figure 4. Without G1 recovery steps, the optimised shape gradually deviates from the initially exact G1 continuity and introduces an average error of around 1.1ı between the tangent plane normals after only 10 design steps. With the G1 recovery steps applied after each cost-function-based perturbation, the deviation remains negligible. 7. OPTIMISATION RESULTS AND DISCUSSION The NsPCC method is applied to the 3D segment of an S-Bend air duct shown in Figure 5. This testcase has been provided by Volkswagen AG to the FlowHead research project (see

Figure 5. Initial geometry consisting of 30 NURBS patches viewed in CATIA V5. Copyright © 2013 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Fluids (2013) DOI: 10.1002/fld

S. XU, W. JAHN AND J.-D. MÜLLER

Acknowledgements). The goal is to deform the central S-section such as to minimise the massaveraged total pressure loss defined as R R E  nE dS C out let Pt ot uE  nE dS i nlet Pt ot u R . (32) JD E  nE dS i nlet u An inlet velocity of 0.1 m/s, a zero back pressure and no-slip walls are used as boundary conditions. The Reynolds number using the inlet height as the reference length is Re D 300. The initial geometry is shown in Figure 5; the duct bends upwards and sideways, so the optimisation result is also expected to be asymmetric with respect to the vertical plane. Only the patches in the cranked S-section of the duct are allowed to move in order to represent build constraints for the inlet and outlet sections. The fixed patches at inlet and outlet are G0 continuous, and G1 continuity is initially satisfied and imposed across interfaces of the moveable patches. The deformable S-section consists of eight NURBS patches in total. Four of those patches each have 96 (16  6) control points and the other four each have 64 (16  4) control points, resulting in a total of 640 control points for the S-section, equivalent to 1920 DoF as each control point is allowed to move in the x, y and ´ directions. At common edges between the eight moveable and fixed NURBS patches upstream and downstream of the S-section, the first three layers of control points of the moveable NURBS patches are fixed so that the entry and exit throats both have G2 continuity and meet with zero curvature. In summary, a total of 400 control points are allowed to move in three directions, equivalent to 1200 DoF. To impose continuity constraints, a total of 240 pairs of test points are distributed along the eight deformable joint edges, 30 on each edge, resulting in 1440 constraint equations. The number 30 is determined as described in Section 5, and this choice leads to a null space with around 570 allowable modes, which may vary slightly over optimisation iterations because of the nonlinearity of G1 constraint function. The corresponding computational surface mesh of the S-section has a total of 3840 boundary nodes. The steepest descent method is used as an optimiser, and the step size is determined using a simple line search bounded by a maximum allowable step size, which requires that the maximum surface node perturbation does not exceed the minimum spacing of the surface mesh. This ensures that the surface mesh perturbation after each design iteration is bounded and the volume mesh smoothing algorithm remains stable. The convergence history of the cost function is shown in Figure 6. As shown in Figure 4, and also confirmed by manually checking the optimised geometry in the CAD program CATIA V5, the final shape remains G1 continuous across patch interfaces where specified. The cost function drops about 25.5% after 60 design optimisation iterations; during each iteration, five normal steps are used to recover G1 continuity. The initial shape and the optimised shape after 60 iterations are shown in Figure 7. The optimised shape shows a distinct bulge at the top and bottom of the S-section, as well as along the four edges, while hollowing inwards from both sides. Also shown in Figure 7 is the surface sensitivity of the

Cost function / Pa

0.016

0.015

0.014 Improvement 25.5%

0.013

0.012

0.011 0

10

20

30

40

50

60

Optimisation iterations

Figure 6. Cost function convergence history. Copyright © 2013 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Fluids (2013) DOI: 10.1002/fld

CAD-BASED SHAPE OPTIMISATION WITH CFD USING A DISCRETE ADJOINT

Initial

Left view

Right view

Optimised

Left view

Right view

Figure 7. Comparison between the initial (upper) and optimised (lower) shapes and the surface sensitivity map for both shapes.

Figure 8. Contour plots of velocity magnitude for the initial (left) and optimised (right) ducts.

Figure 9. Streamlines plots for the initial (left) and optimised (right) ducts. Colour along the streamlines is for illustrative purpose, not related to any flow variable.

initial and the optimised shapes. The sensitivity of the cost function with respect to the normal displacement of the surface is reduced to almost zero everywhere expect at the inlet and outlet throats where the surface is not allowed to move. The resulting shape changes are not intuitive; hence, to better understand the optimisation result, we compare the flowfields. Figures 8 and 9 show the contour plots of velocity magnitude and the streamlines for the initial and optimised shapes. It can be seen that the flowfield of the initial shape has very strong secondary flows, and the cross-sectional cuts exhibit strong non-uniformity in the flow speed. This phenomenon of secondary flow in bent ducts is well known as Dean vortices [34]. There is a large area of separated flow at the bottom of the outlet section in the nearside corner, which will add to the total pressure loss. The flowfield of the optimised shape on the right of Figures 8 and 9 exhibits much weaker secondary flows. The streamlines and velocity magnitude contours show that the separation is significantly reduced. The main difference compared with the original flow lies in the suppressed Dean vortices. Copyright © 2013 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Fluids (2013) DOI: 10.1002/fld

S. XU, W. JAHN AND J.-D. MÜLLER

Figure 10. Cross section view of the S-section for the initial and optimised S-Bend.

Figure 11. Schematics illustrating the formation of strake-like shape. Left: original cross section with a pair of vortices whirling in counter directions; middle: two strakes can be added to suppress the Dean vortices to minimise energy loss; right: alternatively, the sides can hollow inward to mimic the effect of stakes. Table I. Computational cost breakdown for different chains of the optimisation loop. Percentage over total time Primal (Res D 1E5) Adjoint (Res D 1E5) SVD null space G1 recovery steps Shape perturbation

17.28 78.67 2.95 1.01 0.08

It is apparent that this suppression is achieved by the hollowing in of both sides of the S-section, which resembles strakes to suppress the formation of vortices, such as those widely used on aeroplanes, ship hulls and pipelines. The strake-like shape is more apparent in Figure 10, which shows a cross cut in the S-section viewed from the outlet into the S-Bend. Ideally, a pair of stakes along the side of the S-section (Figure 11, middle) are to be formed to suppress the Dean vortices, as is shown on the left of Figure 11. However, as shape optimisation methodology such as NsPCC does not allow a change of topology, the optimised shape (Figure 11, right) clearly shows the tendency to form two strakes on both sides to suppress the secondary flow due to Dean vortices. The whole optimisation process takes a total of 60 optimisation steps on a mesh with 250 000 hexahedra with the dominant cost arising from the primal and adjoint solves. The main steps of the optimisation loop are (a) the primal and adjoint computation and (b) the SVD at each optimisation step to compute the null space of the continuity constraint and thus to find the shape displacement modes, as well as (c) the G1 recovery steps and (d) the shape perturbation step. The costs of each of these steps are shown in Table I. It can be seen that the steps related to the CAD-based parametrisation (b)–(d), that is, SVD null space, G1 recovery steps and the shape perturbation, have a negligible computational cost compared with the primal and adjoint computations. 8. SUMMARY The novel NsPCC shape optimisation methodology has been presented that allows to include the modification of a CAD description inside the design loop. The method uses the NURBS control Copyright © 2013 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Fluids (2013) DOI: 10.1002/fld

CAD-BASED SHAPE OPTIMISATION WITH CFD USING A DISCRETE ADJOINT

points of the generic BRep as design variables and writes the updated set of NURBS control points to a STEP file. The parametrisation is hence fully automatic, not requiring the user to set up a parametrisation in a proprietary CAD system. Seamless connection from and to CAD software packages is carried out through the standardised STEP file format. Geometric continuity at interfaces between NURBS patches is imposed as constraints throughout the optimisation process. The constraints are numerically evaluated at a set of test points, and the number of required test point pairs is determined by evaluating the change in rank of the constraint matrix using SVD. A well-conditioned orthogonal basis for the design space arises from this SVD, allowing to neglect minor modes with little effect on the constraints by appropriately selecting the threshold for singular values considered to be zero. Constraints are maintained linearly by requiring the control point displacements to remain in the null space of the constraint matrix and by perpendicular recovery steps in the range of the constraint matrix to correct for the nonlinear behaviour of G1 and higher constraints. The effectiveness of the method is demonstrated by applying it to the 3D segment of a Volkswagen air duct, achieving a performance improvement of 25.5%. The strake-like shape on the sides of the optimised duct very effectively suppresses the formation of Dean vortices and flow separation and thus effectively reduces the energy loss. The resulting shape shows very strong and detailed deformations with guaranteed smoothness. The proposed method hence satisfies the key requirements on design parametrisations for industrial application: (a) the parametrisation is fully automatic in an open, vendor-neutral form, able to be coupled interchangeably to a range of CAD systems; (b) the design space is as rich as can possibly be expressed in the given BRep and is by construction smooth: the regularity inside patches is given by the order of the patch basis functions, and regularity across patch interfaces is controlled by user-imposed geometric continuity constraints; (c) it permits imposition of geometric constraints such as geometric continuity, and also buildspace, manufacturing or maximum curvature; (d) it returns the optimised geometry in a CAD format for further analysis and processing; and (e) the runtime of the derivative evaluation of the surface node displacements with respect to CAD shape parameters is very low, allowing tight coupling into a simultaneous (one-shot) design loop. ACKNOWLEDGEMENTS

The authors would like to acknowledge the financial support by the European Commission, Grant No. SST.2007-RTD-1 for the FlowHead (Fluid Optimisation Workflows for Highly Effective Automotive Development Processes) project http://flowhead.sems.qmul.ac.uk. REFERENCES 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13.

Pironneau O. On optimum design in fluid mechanics. Journal of Fluid Mechanics 1974; 64:97–110. Jameson A. Aerodynamic design via control theory. Journal of Scientific Computing 1988; 3:233–260. Thévenin D, Janiga G (eds). Optimization and Computational Fluid Dynamics. Springer: Berlin, 2008. Othmer C. A continuous adjoint formulation for the computation of topological and surface sensitivities of ducted flows. International Journal for Numerical Methods in Fluids 2008; 58(8):861–877. DOI: 10.1002/fld.1770. Giannakoglou KC. Design of optimal aerodynamic shapes using stochastic optimization methods and computational intelligence. Progress in Aerospace Sciences 2002; 38(1):43–76. Giles MB, Duta MC, Müller J-D, Pierce NA. Algorithm developments for discrete adjoint methods. American Institute of Aeronautics and Astronautics Journal 2003; 41(2):198–205. Griewank A, Walther A. Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation, second edition. SIAM: Philadelphia, PA, 2008. Giering R, Kaminski T, Slawig T. Generating efficient derivative code with TAF: adjoint and tangent linear Euler flow around an airfoil. Future Generation Computer Systems 2005; 21(8):1345–55. Hascoët L, Pascual V. Tapenade 2.1 user’s guide. Technical Report 0300, INRIA, 2004. Walther A, Griewank A, Vogel O. ADOL-C: automatic differentiation using operator overloading in C++. Proceedings in Applied Mathematics and Mechanics 2004; 2(1):41–44. Cusdin P, Müller J-D. On the performance of discrete adjoint CFD codes using automatic differentiation. International Journal for Numerical Methods in Fluids 2005; 47(6-7):939–945. Jones D, Christakopoulos F, Müller J-D. Preparation and assembly of adjoint CFD codes. Computers and Fluids 2011; 46(1):282–286. Samareh JA. Survey of shape parameterization techniques for high-fidelity multidisciplinary shape optimization. American Institute of Aeronautics and Astronautics Journal 2001; 39(5):877–884.

Copyright © 2013 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Fluids (2013) DOI: 10.1002/fld

S. XU, W. JAHN AND J.-D. MÜLLER

14. Fudge DM, Zingg DW, Haimes R. A CAD-free and a CAD-based geometry control system for aerodynamic shape optimization. American Institute of Aeronautics and Astronautics CP 05-0451 2005. 15. Jameson A, Vassberg A. Studies of alternative numerical optimization methods applied to the brachistochrone problem. Computational Fluid Dynamics Journal 2000; 9(3). 16. Jaworski A, Müller J-D. Toward modular multigrid design optimisation. Lecture Notes in Computational Science and Engineering 2008; 64:281–291. 17. Jakobsson S, Amoignon O. Mesh deformation using radial basis functions for gradient-based aerodynamic shape optimization. Computers and Fluids 2007; 36(6):1119–1136. 18. Reuther J, Jameson A, Alonso JJ, Remlinger MJ. Constrained multipoint aerodynamic shape optimisation using an adjoint formulation and parallel computers, part 1. Journal of Aircraft 1999; 36(1):51–60. 19. Samareh JA. Aerodynamic shape optimization based on free-form deformation. American Institute of Aeronautics and Astronautics Paper; 2004 – 4630:2004. 20. Robinson TT, Armstrong CG, Chua HS, Othmer C, Grahs T. Sensitivity-based optimization of parameterised CAD geometries. 8th World Congress on Structural and Multidisciplinary Optimisation, Lisbon, 2009. 21. Martn MJ, Andrs E, Widhalm M, Bitrin P, Lozano C. Non-uniform rational B-splines-based aerodynamic shape design optimization with the DLR TAU code. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering 2012; 226(10):1225–1242. 22. Schmidt S, Schulz V. Impulse response approximations of discrete shape Hessians with application in CFD. SIAM Journal on Control and Optimization 2009; 48(4):2562–80. 23. Nemec M, Aftosmis MJ. Adjoint algorithm for cad-based shape optimization using a Cartesian method. American Institute of Aeronautics and Astronautics Paper; 2005 – 4987:2005. 24. Robinson TT, Armstrong CG, Chua HS, Othmer C, Grahs T. Optimizing parameterized CAD geometries using sensitivities based on adjoint functions. Computer-Aided Design & Applications 2012; 9(3):253–268. 25. Hazra SB, Schulz V, Brezillon J, Gauger NR. Aerodynamic shape optimization using simultaneous pseudotimestepping. Journal of Computational Physics 2005; 204:46–64. 26. Jaworski A, Cusdin P, Müller J-D. Uniformly converging simultaneous time-stepping methods for optimal design. In Eurogen, et. al. R. Schilling (ed.). Eccomas: Munich, 2005. 27. Yu G, Müller J-D, Jones D. CAD-based shape optimisation using adjoint sensitivities. Computers and Fluids 2011; 46(1):512–16. 28. Piegl L, Tiller W. The NURBS Book, 2nd edition. Springer: New York, NY, 1997. 29. Zhang Q, Lee J-D, Wendisch J-H. An adjoint-based optimization method for helicopter fuselage backdoor geometry. Thirty-Sixth European Rotorcraft Forum, Paris, France, 2010. 30. Ferziger JH, Peri´c M. Computational Methods for Fluid Dynamics. Springer: New York, NY, 2002. 31. van der Vorst H. Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM Journal on Scientific and Statistical Computing 1992; 13(2):631–644. 32. Blom F. Considerations on the spring analogy. International Journal for Numerical Methods in Fluids 2000; 32(6):647–668. 33. Christakopoulos F, Jones D, Müller J-D. Pseudo-timestepping and validation for automatic differentiation derived code. Computers and Fluids 2011; 46(1):174–179. 34. Berger SA, Talbot L, Yao LS. Flow in curved pipes. Annual Review of Fluid Mechanics 1983; 15(1):461–512.

Copyright © 2013 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Fluids (2013) DOI: 10.1002/fld